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^ , Abstract. We study the thermodynamics of quantum particles with long-range interactions at T = 0. 

' Specifically, we generalize the Hamiltonian Mean Field (HMF) model to the case of fermions and bosons. 

' In the case of fermions, we consider the Thomas-Fermi approximation that becomes exact in a proper 

, thermodynamic limit. The equilibrium configurations, described by the Fermi (or waterbag) distribution, 

are equivalent to polytropes with index n — 1/2. In the case of bosons, we consider the Hartree approxi- 
mation that becomes exact in a proper thermodynamic limit. The equilibrium configurations are solutions 
' fH ' , of the mean field Schrodinger equation with a cosine interaction. We show that the homogeneous phase, 

^ ' that is unstable in the classical regime, becomes stable in the quantum regime. This takes place through a 

^ ' first order phase transition for fermions and through a second order phase transition for bosons where the 

S , control parameter is the normalized Planck constant. In the case of fermions, the homogeneous phase is 

stabilized by the Pauli exclusion principle while for bosons the stabilization is due to the Heisenberg un- 
certainty principle. As a result, the thermodynamic limit is different for fermions and bosons. We point out 
analogies between the quantum HMF model and the concepts of fermion and boson stars in astrophysics. 
^ ' Finally, as a by-product of our analysis, we obtain new results concerning the Vlasov dynamical stability of 

the waterbag distribution. We show that spatially homogeneous waterbag distributions are Vlasov stable 

eifF e > Cc = 1/3 and spatially inhomogeneous waterbag distributions are Vlasov stable iff e < e* = 0.379 
' and b > bt = 0.37 where e and b are the normalized energy and magnetization. The magnetization curve 

^ ' displays a first order phase transition at et — 0.352 and the domain of metastability ranges from to e* 

o 



PACS. 5.30.-d Quantum statistical mechanics - 05.45.-a Nonlinear dynamics and chaos - 05.20.Dd Kinetic 



^ I theory - 64.60.De Statistical mechanics of model systems 

^ , 1 introduction classical and quantum self-gravitating systems. This will 

■ place our study in a more general perspective. 

OO . The dynamics and thermodynamics of systems with long- The statistical mechanics of classical self-gravitating 

range interactions (self-gravitating systems, geophysical systems, such as elliptical galaxies and globular clusters, 

, flows, non neutral plasmas,...) has recently received a par- started with the works of Ogorodnikov |llj . Antonov |12j 

PsJ ■ ticular attention from the community of statistical me- and Lynden-Bell & Wood [13] (see reviews |14I15I16] ). 

_ chanics [1121314] . Surprisingly, long-range interacting sys- These authors studied the equilibrium configurations of 

■ terns had not been studied at a general level until now, a self-gravitating gas of classical particles enclosed within 
although the main concepts (such as ensemble inequiv- a box so as to prevent its complete evaporation. They 
alence, negative specific heats, quasi stationary states, considered the microeanonieal ensemble and discovered 

T ' violent coUisionless relaxation and slow collisional relax- the important phenomenon of "gravothcrmal catastro- 

r^J I ation) had been understood early in astrophysics and phe" which takes place below a critical energy Ec = 

^ , two-dimensional turbulence (see, e.g., |5] and references —0.335GM^ / R. This can be viewed as a sort of phase 

therein). Recently, these fundamental concepts have been transition between a gaseous state and a "singular" state 

illustrated and emphasized in the framework of a sim- in which a tightly bound binary star is surrounded by a hot 

pie toy model of systems with long-range interactions halo of stars. We call it "singular" because the maximum 

called the Hamiltonian Mean Field (HMF) model [5]. This entropy state is reached when two stars approach each 

model displays many analogies with self-gravitating sys- other with no limit (leading to infinite density) so that the 

terns |7l8l9ll0j which have proven to be very fruitful. Until potential energy tends to — oo and the temperature of the 

now, only the classical HMF model has been considered halo increases to -l-oo so as to allow for the overall conser- 

In this paper, developing further the analogy with vation of energy. This core-halo state has infinite entropy, 

self-gravitating systems, we shall give a first discussion of On the other hand, Kiessling [T7| and Chavanis [18] con- 

the quantum HMF model for fermions and bosons. Before sidered the canonical ensemble and investigated the phe- 

that, we briefly review the main steps in the history of nomcnon of "isothermal collapse" which takes place be- 
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low a critical temperature ksTc — GMm/2.52R. In that 
case, the final outcome of the collapse is a Dirac peak con- 
taining all the particles. This compact object has infinite 
free energy. These singular states, binaries stars in the 
microcanonical ensemble and Dirac peaks in the canoni- 
cal ensemble, are obtained by assuming that the particles 
are point-like and that quantum mechanics effects can be 
neglected. 

The case of quantum particles in gravitational inter- 
action has also been considered in astrophysics and led to 
the concepts of fermion and boson stars. 

Fermion stars have found astrophysical applications in 
the context of white dwarf stars, neutron stars and dark 
matter models made of massive neutrinos |19I20| . Their 
equilibrium results from a balance between the gravita- 
tional attraction and the quantum pressure due to the 
electrons (in white dwarfs) , the neutrons (in neutron stars) 
or the neutrinos (in dark matter). Gravitational collapse is 
therefore prevented by the Pauli exclusion principle as first 
realized by Fowler [21]. At T = 0, the gas is completely 
degenerate and, in the non relativistic limit, the system 
is equivalent to a poly trope of index n = 3/2. The den- 
sity profile has a finite support leading to a configuration 
with a well-defined radius R. The mass-radius relation of 
classical fermion stars is given 

When special [22] or general [23] relativity is taken into 
account, it is found that no equilibrium state exists above 
a maximum mass, called the Chandrasekhar mass, scal- 
ing like Mch ^ M^lm^ where Mp = {Tlc/GY/'^ is the 
Planck mass. Massive stars cannot pass into the white- 
dwarf stage and become neutron stars. Even more mas- 
sive stars undergo gravitational collapse and become black 
holes [50]. The case of self-gravitating fermions at non- 
zero temperature has been considered by several authors 
|24I25I26I27I^ . The equilibrium states of the non rela- 
tivistic Fermi gas are obtained by coupling the Fermi- 
Dirac statistics to the Poisson equation (this corresponds 
to the Thomas- Fermi approximation). In that case, the 
system extends to infinity and one is forced to enclose the 
particles within a finite box to prevent their evaporation. 
The shape of the caloric curve T{E) depends on a dimen- 
sionless parameter ^ which can be viewed as a normalized 
system size or as an inverse normalized Planck constant 
|16| . For — !■ +CX3, one recovers the classical limit in which 
the system undergoes gravitational collapse and forms a 
singularity. However, when quantum mechanics is taken 
into account, the gravitational collapse stops when the 
system feels the Pauli exclusion principle. In that case, it 
ends up on a nonsingular equilibrium state which typically 
has the form of a completely degenerate and very compact 
nucleus (fermion ball) surrounded by a hot and almost 
homogeneous atmosphere (halo). One can therefore de- 
scribe interesting zeroth and first order phase transitions 
between "gaseous" and "condensed" states, and evidence 
microcanonical and canonical critical points that differ in 
the two ensembles. The complete phase diagram of the 
self-gravitating Fermi gas is given in [16] . 

Boson stars were introduced by Kaup [IH] and Ruffini 
& Bonazzola [^ in the sixties although no astrophysical 



application of these objects was known at that time. Later 
on, it was suggested that dark matter could be made of 
bosons and that boson stars could have formed by Jeans 
gravitational instability [3T]. At T = 0, bosons form a 
Bose-Einstein condensate (BEG) and they are described 
by a single wavefunction ■0(r,t). In the nonrelativistic 
Newtonian limit, the structure of a self-gravitating BEG is 
obtained by solving the Schrodinger-Poisson system and 
in the relativistic limit one must couple the Klein-Gordon 
equation to the Einstein field equations. The Newtonian 
approximation is valid for sufficiently small masses and 
yields the mass-radius relation MR = 9.9h^/Gm? [50] . 
The radius decreases as mass increases, like for classi- 
cal fermion stars, but the scaling is different. This rela- 
tion is valid as long as the radius is much larger than 
the Schwarzschild radius Rs = 2GM/(?. When relativis- 
tic effects are taken into account, there exists a maxi- 
mum mass, the Kaup mass MKaup = 0.63 3M|/m , above 
which no equilibrium configuration exists [29130] . In that 
case, the system collapses to a black hole. Below the 
critical mass, the gravitational collapse of boson stars is 
prevented by the Heisenberg uncertainty principle while 
the gravitational collapse of fermion stars is prevented 
by the Pauli exclusion principle. This is why the Kaup 
mass MKaup ~ Mp/m for boson stars is much smaller 
than the Ghandrasekhar mass Mch ^ Mp jw? for fermion 
stars by a factor m/Mp <C 1. Such small masses (e.g. 
MKaup ~ 10"fcg ~ IO-i^Mq for m - IGe^/c^!) led to 
the belief that boson stars are not very relevant astro- 
physical objects. However, the maximum mass of boson 
stars can be considerably increased by taking into ac- 
count the self-interaction of the particles via a po- 
tential [32] . In that case, the maximum mass scales like 
M ~ X^/'^Mch and becomes comparable with the Chan- 
drasekhar mass of self-gravitating fermions when A ~ 1. 
Therefore, self-coupling can significantly change the phys- 
ical dimensions of boson stars, making them much more 
astrophysically interesting. Recently, some authors [55] 
have proposed that dark matter galactic halos themselves 
could be gigantic self-gravitating Bose-Einstein conden- 
sates (BEG). At the galactic scale, it is possible to neglect 
relativistic effects and use the Newtonian approximation. 
When self-coupling is taken into account, the structure of 
a self-gravitating BEG is obtained by solving the Gross- 
Pitaevskii- Poisson system (this corresponds to the Hartree 
approximation). The short-range interaction produces an 
effective pressure described by a barotropic equation of 
state equivalent to a polytrope of index n = 1. The equi- 
librium of the system then results from the balance be- 
tween the gravitational attraction, the pressure due to 
the self-interaction and the quantum pressure due to the 
Heisenberg uncertainty principle. A detailed study of the 
equilibrium configurations is made in [34| . 

In this paper, we shall apply the same concepts to the 
HMF model in which the gravitational interaction in d = 3 
dimensions is replaced by a cosine interaction in d = 1 di- 
mension. This generalization has two main interests. First, 
it will allow us to study the thermodynamics of quantum 
particles with long-range interactions in a simpler model 
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and to derive several analytical results. Secondly, it will be 
possible in future works to perform very efficient numeri- 
cal simulations of this system since the HMF model is one 
dimensional and the potential of interaction is smooth. 
This will allow us to explore the dynamical properties of 
quantum particles with long-range interactions in a sim- 
ple setting. Generalizing the HMF model to the quantum 
regime is the natural next step in the systematic explo- 
ration of the properties of this model since its introduction 
in 1995 [6]. 

The statistical mechanics of the classical HMF model 
has been treated by different methods (see review [3]). 
Among them, the approach of |7I9I35| based on the maxi- 
mization of the Boltzmann entropy at fixed mass and en- 
ergy closely follows the standard methodology developed 
in astrophysics |12I13I14I15ITB] . For the classical HMF 
model, there exists a second order phase transition be- 
tween a homogeneous phase (for E > Ec and T > Tc) and 
a clustered phase (for E < Ec and T < Tc). Contrary to 
3D self-gravitating systems, the ensembles are equivalent 
for all accessible energies and temperatures. At T = or 
E = Emim a limit that will be considered in the following, 
the density profile p{6) = M6{6) is a Dirac peak centered 
ai 9 = (say) containing all the mass. This "singular" 
state corresponds to a complete condensation of the sys- 
tem in which the maximum magnetization is 6 = 1. We 
shall here study how this result is modified when quantum 
mechanics effects are taken into account. 

In the case of fermions, the statistical equilibrium state 
is obtained by maximizing the Fermi-Dirac entropy at 
fixed mass and energy. This leads to the mean field Fermi- 
Dirac distribution with a cosine interaction. In fact, this 
problem has already been treated in a different context. 
Indeed, in the coUisionless regime of the dynamics, the 
classical HMF model is governed by the Vlasov equation. 
Now, the Vlasov equation can undergo a process of vio- 
lent coUisionless relaxation leading to a quasi stationary 
state (QSS). In the case where the fine-grained distribu- 
tion function takes only two values f = rjo and / = 0, 
the statistical mechanics of the Vlasov equation developed 
by Lynden-Bell [33] predicts that the QSS is obtained by 
maximizing a "fermionic" mixing entropy at fixed mass 
and energy. Here, the "degeneracy" is due to dynami- 
cal effects (Liouville theorem), not to quantum mechanics 
(Pauli's exclusion principle). Still, the mathematical prob- 
lem is the same provided that the maximum value of the 
distribution function /o fixed by the initial condition is in- 
terpreted as the maximum value of the distribution func- 
tion 2/h fixed by the Pauli exclusion principle (see [37] 
in the astrophysical context). The phase transitions asso- 
ciated with the Lynden-Bell or Fermi-Dirac distributions 
are very rich and subtle and they have been described in 
detail in |38l39l40l41l42l43l44j . Due to the analogy with 
the quantum problem, these studies directly apply to a gas 
of fermions with cosine interaction. Therefore, we shall not 
repeat this analysis here. Wc shall, however, study the case 
T = (corresponding to the ground state energy Eground) 
that was not treated in detail in the previous works. In 
that limit, the system is completely degenerate and the 



distribution function reduces to the Fermi distribution or, 
equivalently, to the (spatially inhomogeneous) waterbag 
distribution. This distribution is equivalent to a polytrope 
of index n = 1/2. The general theory of polytropes in the 
context of the HMF model has been recently developed 
by Chavanis & Campa [3S] and we shall make use of their 
results in the specific case n = 1/2 which presents interest- 
ing features. In terms of the normalized Planck constant x 
for fermions defined in equation (|77p , we show that the ho- 
mogeneous phase is unstable for x < Xc = V^, metastable 
for Xc < X < Xt — 1-45 and fully stable for x > Xt- In Par- 
allel, the inhomogeneous phase is fully stable for x < Xt, 
metastable for xt < X < X* — 1-48 and disappears for 
X > X* (see Figure (TS]). This shows that quantum ef- 
fects can stabilize the homogeneous phase because of the 
Pauli exclusion principle, exactly like for fermion stars in 
astrophysics. For the HMF model, this stabilization hap- 
pens through a first order phase transition where the con- 
trol parameter is the normalized Planck constant. As a by 
product of our analysis, we obtain new results concerning 
the Vlasov dynamical stability of the waterbag distribu- 
tion. We show that spatially homogeneous waterbag dis- 
tributions are Vlasov stable iff e > Ec = 1/3 and spatially 
inhomogeneous waterbag distributions are Vlasov stable 
iff e < = 0.379 and > fe, = 0.37 where e and b are the 
normalized energy and magnetization. The magnetization 
curve displays a first order phase transition at = 0.352 
and the domain of metastability ranges from ec to e* (see 
Figure g]). 

At T = 0, bosons form a Bose-Einstein condensate 
(BEC) and the equilibrium state is obtained by solving the 
mean field Schrodinger equation with a cosine interaction 
(the pressure derived from the Bose-Einstein statistics is 
zero at T = so that a hydrodynamical description is 
not appropriate). The normalized Planck constant x ioi 
bosons is defined in equation (|126[) . For x = 0, we re- 
cover the classical result: the density profile forms a Dirac 
peak at 9 = 0. For x — ^ 0, we can make a harmonic ap- 
proximation and obtain quantum corrections at the order 
0(x)- More generally, by solving the problem numerically 
for any value of x, we evidence a second order phase tran- 
sition between magnetized states for X < Xc = and 
non magnetized states for x > Xc- This shows that quan- 
tum effects can stabilize the homogeneous phase because 
of the Hcisenberg uncertainty principle, exactly like for 
boson stars in astrophysics. 

The paper is organized as follows. In Section [51 we 
introduce the HMF model. In Section [3] we consider the 
case of fermions at T = and study the cosine Fermi 
distribution. We show that our study also determines 
the Vlasov dynamical stability of the (possibly inhomo- 
geneous) Lynden-Bell and waterbag distributions, inde- 
pendently on the quantum mechanics context. In Section 
m we consider the case of bosons at T = and study 
the cosine Schrodinger equation. In Section [5] we con- 
sider the mean field Gross-Pitaevskii equation describing 
BECs with short-range interactions or fermions beyond 
the Thomas-Fermi approximation. In this context, we 
study the stability of a homogeneous distribution with re- 
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spect to the quantum Eulcr and Vlasov (or Wigncr) equa- 
tions. The appropriate thermodynamic hmits for fermions 
and bosons are discussed in Appendix]^ Technical results 
are given in the other Appendices. 



2 The HMF model 

Wc consider a system of N particles of unit mass m = 1 
moving on a ring of unit radius R = I and interacting via 
a cosine potential of the form u = — ^ cos{0 — 9'), where 
fc > is the coupling constant. The Hamiltonian reads 



N 



H 



1 ^ — A 2 ^ \ — ^ 



i=i 



(1) 



where 9i and Vi — 6i are the position (angle) and velocity 
of particle i. We introduce the magnetization vector b = 
{bx,by) where = j^J^i cos 9i and &y = ;^ I]i sinflj. 

We assume that the system can be described by a dis- 
tribution function f{9,v,t) such that f{9,v,t)d9dv gives 
the density of particles with position 9 and velocity v at 
time t. It is normalized such that AI — J f d9dv. As we 
shall see, this description applies to classical particles and 
fermions (in the Thomas-Fermi approximation) but not 
to bosons at r = forming Bose-Einstein condensates. 
In the mean field approximation, the energy (kinetic -|- 
potential) is given by 

E = K + W = f d9dv + ^ / P^d9, (2) 



where 



k 



2ir 



cos{9 ~9')p{9\t)d9', 



(3) 



is the self-consistent potential generated by the density 
of particles p{9,t) = J f{9,v,t)dv. Expanding the cosine 
function, the potential can be rewritten 



<P{9, t) = -Bx cos 9- By sin 9, 



where 



Bx 



B,, 



2tt 



2^ JO 



p{9,t) cos 9 d9, 



pi9,t) sm9d9, 



(4) 

(5) 
(6) 



are proportional to the two components of the average 
magnetization bx = J p cos 9 d9 and by = j psm9d9. 
The potential energy can be expressed in terms of the 
magnetization as 

On the other hand, the kinetic energy can be expressed in 
term of the pressure p{9, t) = J fv^ dv as 



E : 



pd9. 



(8) 



3 Fermions at zero temperature 

In this section, we consider the HMF model at T = in 
the case where the particles are fermions of spin s = 1/2. 
We use a mean field approximation that becomes exact 
in a proper thermodynamic limit N — > +oo defined in 
Appendix |21 This mean field approximation could be rig- 
orously justified like in the astrophysical problem |24l25j . 



3.1 The Fermi-Dirac distribution 

According to the Pauli exclusion principle, there are at 
most 2s -|- 1 = 2 fermions in a phase space cell of size 
h, where h is the Planck constaniQ. Therefore, the maxi- 
mum value of the distribution function fixed by the Pauli 
exclusion principle is 



2 
h' 



(9) 



The statistical equilibrium state in the microcanonical en- 
semble is obtained by maximizing the Fermi-Dirac entropy 



S = - 



Vo Vo 



1 - ^ 1 In 



Vo 



d9dv, 



(10) 

at fixed mass and energy (see, e.g., [16]). Writing the vari- 
ational principle in the form 5S—/35E—a5M = 0, where /3 
and a are Lagrange multipliers, we obtain the mean field 
Fermi-Dirac distribution 



1 + Ae'3''o(T+*(»)) 



(11) 



This distribution function satisfies the constraint / < '/o- 
The statistical equilibrium state in the canonical ensem- 
ble is obtained by minimizing the Fermi-Dirac free energy 
F ~ E — TS at fixed mass, where T > is the tem- 
perature. The critical points of constrained entropy and 
constrained free energy are the same, given by equation 
(jlip . but their stability may differ in microcanonical and 
canonical ensembles in case of ensemble inequivalence (see, 

e.g., m)- 



3.2 Tlie Fermi distribution 



At T = 0, the Fermi gas is completely degenerate since 
the states with individual energy e = w^/2 -t- ^{9) smaller 
than the Fermi energy ep are completely filled. This leads 
to the Fermi distribution 



f = Vo, if e < CF, 
/ = 0, if e>eF. 



(12) 



^ Since the units of length and mass have been fixed to unity, 
the Planck constant h so defined is dimensionless (see Ap- 
pendix |X] for the restoration of dimensional variables). There- 
fore, it can be viewed as an external parameter that can take 
values between (classical regime) and -f cxj (quantum regime). 
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It can be rewritten 



3.3 The homogeneous polytrope n =1/2 



/ = %, if \v\<VFi0), 
/ = 0, if \v\>vf{9), 



where 



(13) 



(14) 



is the space dependent Fermi velocity. This zero temper- 
ature hmit corresponds to the ground state of the Fermi 
gas (minimum energy state Emm)- 

The density and the pressure corresponding to the dis- 
tribution function (jl3|) are given by 



pie) = \ fiv)dv^2i^oVFiO), 

J-VF(d) 



(15) 



fivydv = -rjov%i9). (16) 
-vp{e) 'i 

Ehminating the Fermi velocity between these two expres- 
sions, we obtain the equation of state 



1 



12%^' 



This is the equation of state of a polytrope 



p = Kp^, 7=1 + -, 
n 



with a polytropic constant 



12??^ 12 



and a polytropic index 



7 = 3, i.e. n = -. 

2 



(17) 



(18) 



(19) 



(20) 



As is well-known in astrophysics, the Fermi gas at T = 
(e.g. a white dwarf star) is equivalent to a polytrope with 
index n = d/2 (in our case, the dimension of space is d = 
1) [in]- Accordingly, wc can study the Fermi gas at T = 
by using the theory of polytropes developed by Chavanis 
& Campa [45] in the context of the HMF model. In this 
analogy, the polytropic constant K can be interpreted as 
a polytropic temperature. 

Remark: the Fermi distribution p2[) is a particular 
steady state of the Vlasov equation called the waterbag dis- 
tribution. Therefore, as a by-product, our study will also 
determine the structure and the stability of the (spatially 
inhomogeneous) waterbag distribution, independently of 
the quantum mechanics interpretation. 



Let us first consider the case of a spatially homogeneous 
distribution. In that case, = 0, and the Fermi distribu- 
tion can be rewritten 



/ = m, if kl < VF, 
/ = 0, if \v\>VF, 



(21) 



with vf = \f2tF- The density is given by p = 2riQVF- The 
Fermi velocity is related to the mass and to the maximum 
value of the distribution function by 



Vf 



M 



Mh 



The energy is 



-4 



M 

^ * 27r 



Using equation ((T^ . wc obtain 
E 



967r2772 



96 



(22) 



(23) 



(24) 



If we define the normalized energy and the normalized 
inverse polytropic temperature by |45j : 



^TlE 



the relation 



can be rewritten 

_ 1 



(25) 



(26) 



Of course, the magnetization vanishes in the homogeneous 
phase: & = 0. 



3.4 The inhomogeneous polytrope n = 1/2 

Let us now consider the case of spatially inhomogeneous 
distributions. Combining equations (fTS)) and ([H)) . the den- 
sity profile is given by 



p{9) = 27]oV2(ef -<?(^))- 



(27) 



We need to consider two cases, depending on whether eF 
is positive or negative. 



3.4.1 The case e_F > 

Let us first assume that e^? > 0. In that case, defining 
A = 2r\Q^2tF and using equation the density profile 
(l27l) can be rewritten 



p{d) = A 



1/2 



(28) 
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Wc can assume, without loss of generality, that the distri- 
bution is symmetric with respect to the axis = 0. This 
implies that By = and = B. Then, using equation 
(U), we obtain 



with 



p{0) A ( 1 + -a;cos( 



B 



1/2 



(29) 



(30) 



The amplitude A is determined by the mass M = j pdO 
according to 



A = 



M 



27r/3,o(a;) ' 

where we have introduced the integrals 



(31) 



2tt 



1 + 



7- 



1 \~ 
-xcos9 I cos{m9) d9, 



(32) 



that can be interpreted as deformed Bessel functions [45] . 
By symmetry, we can restrict ourselves to the interval 
< 9 < IT. We need to distinguish two cases [IS]. If 
X < Xc = 3/2, the polytrope is incomplete in the sense that 
the density is strictly positive a.t 9 = n. If x > Xc = 3/2, 
the polytrope is complete in the sense that the density 
vanishes at 9 ~ 9r < tt where 



arccos — 



2x 



(33) 



In that case, the profile has a compact support: p = for 
9c < 9 < TT. Some typical density profiles are represented 
in Figure [T] Since vp{9) = 27^p{9), the density profile 
also gives the shape of the Fermi distribution in the phase 
space {9,v). The velocity distribution (l>{v) = J f d9 is 



2 2 



1)] if-<2|7(l- 



)i/2 and 



(j){v) = 2770 arccos[^(^'' 
0(w) = otherwise. 

Substituting equation (j29p in equation ([5]), we find that 
the magnetization is related to x by 



27tB ^ h.ijx) 
kM hfii^)' 



(34) 



Combining equations (|5n|) , ([5T|) and , we find that the 
normalized inverse polytropic temperature is related to x 

by 



V 



kn 



KM 2hfl{x)h^i{x) 



(35) 



Finally, it is shown in [35] that the normalized energy can 
be written 



SttE ^ 2 /3,1(3:)' _ 4 /3,1(3;) 
kKP 3/3,o(.t)2 3/3,0(2;) 
, 2/3,i(x) 



cos 9r 



X /a, 0(2;) 



(36) 
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Fig. 1. Density profile of n = 1/2 polytropes for different 
values of x. For e_F > we have taken x = 1 and a; = 3. For 
ep < we have taken s = 3. 



For incomplete polytropes [x < Xc), it takes the form 

2 /3.1 (a;)' ^ 4 /3.i(x) ^ 2 /3,i(x) / _ 2^ 
3/3,o(a:)2 3/3,0(2;) a::/3,o(a;) V 3' 



For complete polytropes (a; > Xc), it reduces to 

2 /3,i(x)^ ^ 2/3,1 (x) 
3/3,o(a;)2 xlsfiix)' 



(37) 



(38) 



3.4.2 The case < 

Wc now assume that ep < 0. In that case, defining 
A = 27^0 v2kH ^^'^ using equation (fT9)) . the density pro- 
file ((?7]) can be rewritten 



p{9) = A 



-1 - 



1/2 



'SA^K 

Using equations (|4]) and (pO| . we obtain 

2 



piO) 



-1 + —x cos 9 



1/2 



(39) 



(40) 



The amplitude A is determined by the mass M = J pd9 
according to 



A = 



M 



27rl3,o(x) ' 
where we have introduced the integrals 



(41) 



2ir 



7 - 1 \ 

-IH a; cos 6* I cos{m9)d9. 



(42) 

The central density is defined iff .t > Xc = 3/2, so that 
we shall restrict ourselves to this range of parameters. In 
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that case, the polytropc is always complete in the sense 
that the density vanishes at 6 = 6c < where 

3 



= arecos 



2x 



(43) 



The density has a eompaet support: p = for 6c < < n 

2 2 

otherwise. 



(see Figurc[T]). Furthermore, (j){v) = 27]o arccos[^(^p- 
1)] if^'<4;(-l + HV2and # 

Substituting equation (j40p in equation ([S]), we find that 
the normalized magnetization is related to x by 



27rB 2:3,1(3;) 



kM 



(44) 



2:3,0(2;) ' 

Combining equations (PU)) . (|^T|) and (|33]), the normalized 
inverse polytropic temperature is related to x by 

"^^KM" 2l3,o(a:)2:3,i(a;)' ^ ^ 

Finally, extending the results of [45] , we find that the nor- 
malized energy can be written 



4 2:3,1 (x) 



2 J3,i(a;)^ 

3l3,o(a;)2 3l3,o(x) 

2 13,1 (a;) A 2 

' i X 

3 



a; 2:3,0(2;) V " / + 
Since the poly tropes are always complete, it reduces to 

2X3,1 (x) 



(46) 



2 2:3,i(x)2 



3X3,o(x)2 xl3fi{x)' ^^'^^ 

Remark: the case ej? < was forgotten in |45j where 
it was implicitly assumed that A > in equation (33) 
of that paper. Therefore, it was found that the series of 
equilibria suddenly stops at some finite magnetization, in- 
verse temperature and energy corresponding to x — -|~oo 
in equations (p4)) - ([36| . In fact, if we take the case A < 
into account, the series of equilibria continues until the 
point 6 — 1, 77 — -|-oo and e — >■ e,„i„ = — 2 corresponding 
to x — > Xc = 3/2 in equations (|l^ - (|46p . In that limit, 
the density profile tends to a Dirac peak p{9) = MS{9) 
(see Appendix ID|) . This is more satisfactory on a physical 
point of view. Fortunately, this part of the branch does 
not change the nature of the phase transitions described 
in [45] so that the conclusions of this study are unaltered. 



3.5 Vlasov dynamical stability of the waterbag 
distribution 

Before investigating the thermodynamical stability of the 
Fermi distribution (see Section 13. 6p , we shall make a di- 
gression and investigate its dynamical stability with re- 
spect to the Vlasov equation. In that context, the Fermi 
distribution will be referred to as the waterbag distribu- 
tion. On general grounds, a thermodynamically stable dis- 
tribution must be dynamically stable. Therefore, it makes 
sense to first study its dynamical stability. Furthermore, 
the Vlasov dynamical stability of the waterbag distribu- 
tion is interesting in its own right, independently of the 
quantum mechanics interpretation. 



3.5.1 Formal nonlinear stability 

Let us first recall general results of formal nonlinear dy- 
namical stability. A distribution function of the form 
/ = /(e), depending only on the individual energy e = 
/2 + (l>{9) of the particles, is a steady state of the Vlasov 
equation. If / = /(e) with /'(e) < 0, then it extrem- 
izes a functional S ^ — J C{f) d9dv, where C is convex 
{C" > 0), at fixed mass M and energy E. Furthermore, 
it can be shown that if / is a maximum of S at fixed 
mass and energy, then it is formally nonlincarly dynam- 
ically stable with respect to the Vlasov equation [17]. It 
can also be shown that formal stability implies linear sta- 
bility although the converse is wrong in general. We are 
thus led to investigating the maximization problem 



max {5[/] !£;[/] 



E, AI[f]=M}, 



(48) 



which is similar to a criterion of microcanonical stability 
in thermodynamics. A less refined condition of formal non- 
linear stability is provided by the minimization problem 



nnn{F[f]^E[f]-TS[f] 



I M[f] = M} , (49) 



which is similar to a criterion of canonical stability in ther- 
modynamics. This corresponds to the classical Casimir- 
energy stability criterion |48j . The variational principles 
(|48)) and (|49)) have the same critical points. Furthermore, 
a solution of (|49p is always a solution of the more con- 
strained problem (|15|) but the converse is wrong in case of 
"ensemble inequivalence" |49j (a notion applied here in a 
dynamical context). 

A distribution function of the form / = /(e) with 
/'(e) < determines a "gas" characterized by a density 
profile p{9) and a barotropic equation of state p = p{p) 
[3]. Furthermore, it can be shown |9I50I51) that the min- 
imization problem (|49p is equivalent to the minimization 
problem 



min{F[p] 
p 



M[p] = M} , 



where 



d9 



Pl''P^dp'd9. 



(50) 



(51) 



It can also be shown that this minimization problem pro- 
vides a necessary and sufficient condition of formal non- 
linear dynamical stability of the "gas" with respect to 
the Euler equation |52| . Therefore, a distribution function 
/(e) with /'(e) < is stable with respect to the Vlasov 
equation if the corresponding barotropic gas p{9) is sta- 
ble with respect to the Euler equation, but the reciprocal 
is wrong in case of ensemble inequivalence (since (|49p is 
not equivalent to (j48p '). This provides a new interpretation 
of the (nonlinear) Antonov first law in terms of ensemble 
inequivalence [50] . 

If we consider a functional of the form 



1 



9-1 



ir - I) dOdv, 



(52) 
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then the critical points of "cntropy"Q S at fixed mass and 
energy, satisfying 6S — f35E — aSM = 0, where /3 and 
a are Lagrange multipliers, are given by the polytropic 
distribution function 



q 



-+m 



1/(9-1) 



(53) 



It is convenient to introduce the index n by the definition 

n4 + ^. (54) 

Then, it can be shown [l5] that the equation of state asso- 
ciated with the distribution function (|53p is the polytropic 
one 

Kp-', 7=1 + -, (55) 



P 



where the polytropic temperature K is an increasing func- 
tion of /3~^. On the other hand, the density profile is of 
the form 



A 



7-1 



m 



(56) 



where A can be related to ii. Finally, the free energy ([51 
becomes 



{p' -p)de. 



(57) 



In the context of the HMF model, the Vlasov dynam- 
ical stability of polytropic distribution functions with ar- 
bitrary index n has been studied in |45| . The waterbag 
distribution corresponds to a polytropic index n = 1/2, 
i.e. 7 = 3 and q = oo. Indeed, equation (j53p reduces to 
the waterbag distribution (fT^ . equation ([55| reduces to 
the density profile ([25)1 and equation ([55]) reduces to the 
equation of state ([T7|) . However, the formal stability cri- 
terion based on the maximization problems (|48| or (j49[) 
is not directly applicable since the waterbag distribution 
(fT^ is a singular case for which no associated functional 
S exists (for q = oo the functional (j52p is undetermined). 
In order to investigate its Vlasov dynamical stability, we 
shall use a "ruse" and view the waterbag distribution (|12D 
as a limit of a polytropic distribution with n — > 1/2. Thus, 
we shall use stability results established for polytropes of 
index n > 1/2 and pass to the limit n — 1/2. 

Remark 1: although S[f] and F[f] are ill-defined for 
n = 1/2, the functional F[p\ is well-defined and takes the 
form 



1 



F[p\ ^t; P'^dr 



K 



0^ d9. 



(58) 



We could thus study the Vlasov dynamical stability of 
the Fermi distribution through the minimization problem 
(ISDl). This will be considered in Section [331 where another 



^ As explained in [JS], we use a thermodynamical analogy 
to investigate a dynamical stability problem. By an abuse of 
language we adopt a common vocabulary. 



interpretation of this functional will be given in relation to 
the thermodynamical stability of the Fermi distribution. 

Remark 2: when we study the Vlasov stability of a 
polytrope of index n and mass M, the proper control pa- 
rameter is the energy E (for the maximization problem 
([15[) ) or the polytropic temperature K (for the minimiza- 
tion problem ([49[) ). We shall therefore use these control 
parameters in a first step and express the conditions of 
stability in terms of these parameters. Then, when we fo- 
cus on the waterbag distribution n = 1/2, it is more con- 
venient to take rjo (which is a function of E or K) as a 
control parameter and express the conditions of stability 
in terms of 770. This is particularly important in relation 
to the Lynden-Bell theory where 770 is the natural control 
parameter [38143144] . Finally, when we come back to the 
initial quantum problem, it is more relevant to take h as 
a control parameter (which is related to 770) and express 
the conditions of stability in terms of h. Although these 
problems are closely related, it is important, for clarity, 
to consider them successively. This is what we shall do in 
the sequel. 



3.5.2 Homogeneous distributions 

Let us first consider the Vlasov dynamical stability of 
a spatially homogeneous distribution. We assume that 
/ = f{v) is an even distribution function with a single 
maximum at 1; = (this corresponds to / = /(e) with 
/'(e) < and ^ = 0). Such a distribution is stable iff 



1 



or, equivalent ly, iff 



+00 



fiv) 



>0, 



kM 
1^' 



(59) 



(60) 



where = p'{p) is the velocity of sound in the corre- 
sponding barotropic gas [3]. Equations ([5^ and ([5(I[) are 
equivalent and they provide criteria of linear and nonlin- 
ear dynamical stability with respect to the Vlasov and 
Euler equations j9l47l54l53| . 

For the waterbag (or Fermi) distribution character- 
ized by the equation of state p7]) . we find that Cs{0) = 
p{9)/{2r]o). We note that the velocity of sound coincides 
with the Fermi velocity: Cs{9) = vf{9). For a spatially 
homogeneous distribution p = A//(27r), we obtain 



/ M 



_ flMhV _ 3KAP 



(61) 



Using equations ([6T|) and ([26]) . the stability criterion ([60] 
can be rewritten in terms of the normalized inverse poly- 
tropic temperature and normalized energy as 



V <i]c^ 3, 



e > €c 



(62) 



This is a particular case of the general result (156) of [3] 
corresponding to a polytropic index 7 = 3. Of course, the 
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same results can be obtained directly from the criterion 
((59|) using the distribution function ((2T|) for which f'{v) = 
r]o[S{v + vf) - S{v - vf)]- 



3.5.3 Inhomogeneous distributions 

The Vlasov dynamical stability of spatially inhomoge- 
neous distributions is more difficult to investigate. Sta- 
bility criteria can be obtained by solving eigenvalue equa- 
tions |9I45I53) . variational principles [35147] or dispersion 
relations [SH]- On the other hand, as explained in [5S], 
we can study the formal stability of a steady state of the 
Vlasov equation by using the Poincarc theory of linear se- 
ries of equilibria (see, e.g., |15ll6j ). This is a very simple 
and powerful graphical method that only requires to solve 
the first order variational problem and to know the sta- 
bility of at least one point in the scries of equilibria by 
another method. Let us apply it to the present situation 
(for more details see [45]). 

The formal stability of a steady state of the Vlasov 
equation is based on the optimization problems ()48[) and 
(|49)) . To solve the maximization problem (|48|) correspond- 
ing to the "microcanonical" ensemble, we just need to plot 
/3 (conjugate variable) as a function of E (conserved quan- 
tity) since (3 = {dS/dE)M- It is shown in [45] that 77 is a 
monotonic function of j3 so that it is equivalent to plot 77 
as a function of e. This can be done by eliminating x be- 
tween equations p5p and and between equations (|45p 
and (|46p . The scries of equilibria is represented in Figures 
[2] and [3] The homogeneous phase exists for e > and the 
inhomogeneous phase for e„„;„ = —2 < e < e* ~ 0.379. It 
bifurcates from the homogeneous phase at ec = 1/3. Ac- 
cording to the Poincare theorem, a change of stability can 
occur only at a bifurcation point or at a turning point. 
We have already shown, by a direct calculation, that the 
homogeneous phase is unstable for e < Cc and stable for 
e > Ec (see Section [3. 5. 2p . Therefore, the inhomogeneous 
branch that appears dX e = tc is necessarily unstable. At 
e = there is a turning point of energy so that the 
inhomogeneous branch becomes stable. Since there is no 
other turning point or bifurcation point in the scries of 
equilibria, it remains stable until the end. Wc can use the 
same method in the "canonical" ensemble. To solve the 
minimization problem (|49p . wc just need to plot E as a, 
function of (3 since E — {d{f3F)/dj3)M- As explained pre- 
viously, this is equivalent to plotting e as a function of 
rj. Therefore, we just need to rotate the curves of Figures 
[5] and 13] by 90 degrees. The homogeneous phase exists for 
?7 > and the inhomogeneous phase for 77 > 77* — 2.737. It 
bifurcates from the homogeneous phase at ?7c = 3. We have 
already shown that the homogeneous phase is unstable for 
77 > 77c and stable for rj < rjc (see Section [3.5.2p . There- 
fore, the inhomogeneous branch that appears at t] = rjc is 
necessarily unstable. At 77 = 77*, there is a turning point of 
polytropic temperature so that the inhomogeneous branch 
becomes stable. Since there is no other turning point or bi- 
furcation point in the series of equilibria, it remains stable 



until the end. Wc note that the ensembles are equivalenl|3 
for the polytropic index n = 1 /2 although there is a small 
region of ensemble inequivalence for 71 > 1/2 [5S]. This 
is because, for n = 1/2, the series of equilibria makes a 
spike at (e, 77) = (e,,7;H.) while for n > 1/2 the turning 
points of energy and polytropic temperature are distinct. 
The case tt, = 1/2 is therefore very singular in this respect. 
We shall see below another consequence of this "singular" 
behavior. 




Fig. 2. Series of equilibria of the waterbag distribution corre- 
sponding to a polytrope n = 1/2. It gives the inverse polytropic 
temperature as a function of the energy. The homogeneous 
branch (unstable) tends towards e = and the inhomogeneous 
branch (stable) tends towards e = —2. 
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Fig. 3. Zoom of Figure [5] near the bifurcation point. 



Wc can also study the stability of a polytrope by plot- 
ting the "entropy" S versus the energy E, and by compar- 
ing the entropies of the solutions having the same energy. 
It is shown in [3S| that, for n > 1/2, the normalized en- 

Since the "ensembles" are equivalent, the stability crite- 
ria that we obtain determine both the Vlasov stability of the 
waterbag distribution (|12p and the Euler stability of the cor- 
responding "barotropic gas" described by the density profile 
1271). 
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tropy is given by 

s = - - efcm?7^^, (63) 

where ekin is the normalized kinetic energy. For n > 1/2, 
the study of the curve s(e) is equivalent to the study of 
the curve cr(e) where 

^ = -4-rS'"- (64) 

If we pass to the limit ?i — s- 1/2, we find that 

cr = —rj. (65) 

Therefore, for the polytrope n = 1/2, the entropy a is 
equal to the opposite of the inverse polytropic tempera- 
ture rj. This is a very singular situation. Indeed, the caloric 
curve r]{e) coincides with the entropic curve (T(e)! This is 
the reason why the caloric curve of Figure [31 which has a 
"triangular" shape and presents "cusps" , looks similar to 
the familiar entropy vs energy curve characteristic of first 
order phase transitions (see, e.g.. Figure 9 of [1S|)- There- 
fore, the intersection between the homogeneous and inho- 
mogencous branches in Figure |3] determines the transition 
energy et at which these phases have the same "entropy" . 
Similar results are obtained in the canonical ensemble. In 
that case, we can study the stability of a polytrope by 
plotting the free energy as a function of the polytropic 
temperature K, and by comparing the free energies of the 
solutions having the same polytropic temperature. Now, 
using the results of l45|, we see that for n = 1/2 the free 
energy coincides with the energy, i.e. 

/ = e. (66) 

Therefore, the free energy curve f{r]) coincides with the 
caloric curve e(ry). Therefore, the intersection between the 
homogeneous and inhomogeneous branches in Figure [3] 
determines the transition temperature rjt at which these 
phases have the same "free energy"|3. 

We can now describe more precisely the dynamical sta- 
bility of the watcrbag distribution with respect to the vari- 
ational problems (|^5)) and in the limit n — > 1/2"*". 
In the microcanonical ensemble, a global entropy maxi- 
mum will be called fully stable (S), a local entropy max- 
imum will be called mctastable (M) and a saddle point 
will be called unstable (U). Considering Figure [3l we 
conclude that the homogeneous phase is fully stable for 
e > et — 0.352, metastable for ec < e < e* and unsta- 
ble for < e < Ec. The upper part of the inhomogeneous 
branch is unstable. The lower part of the inhomogeneous 
branch is metastable for and fully stable for 
CmiTi <£<£(. The energy e* can be interpreted as a spin- 
odal point at which the metastable inhomogeneous branch 

* We recall that we are studying here the Vlasov dynamical 
stability problem. Therefore, when we use the word "entropy", 
"free energy", "microcanonical ensemble", "canonical ensem- 
ble" , "temperature" etc. we are invoking the thermodynamical 
analogy explained in [45] ■ To be more correct, we should add 
the prefix "pseudo" in front of these terms. 




Fig. 4. Magnetization as a function of the energy. 




2.5 3 3.5 4 4.5 5 



Fig. 5. Magnetization as a function of the inverse polytropic 
temperature. 



disappears. In the canonical ensemble, a global free energy 
minimum will be called fully stable (S), a local free energy 
minimum will be called metastable (M) and a saddle point 
will be called unstable (U). Considering Figure [3l we con- 
clude that the homogeneous phase is fully stable (S) for 
?/ < ?7t ~ 2.844, metastable (M) for r]t < ?] < rjc and 
unstable (U) for t] > t]c- The upper part of the inhomo- 
geneous branch is unstable. The lower part is metastable 
for rj^ < r] < rjt and fully stable for r] > rjt. The inverse 
temperature 77* can be interpreted as a spinodal point at 
which the metastable inhomogeneous branch disappears. 
As noted previously, the ensembles are equivalent. If we 
regard the curve of Figure[3]as a caloric curve ?7(e), we see 
that r] is continuous at the transition e = while r]'(e) is 
discontinuous. This looks like a second order phase tran- 
sition. However, if we regard the curve of Figure [3| as an 
entropic curve cr(e), we see that it displays a first order 
phase transition marked by the discontinuity of cr'(e) at 
ct and the existence of metastable states. For n > 1/2, 
this leads to a discontinuity of 77(e) at et (see Figure 21 of 
|1S|). Therefore, it is better to say that the caloric curve 
r]{e) displays a first order phase transition (as for n > 1/2) 
in which the jump of temperature tends to zero! The poly- 
trope n = 1/2 (water bag distribution) presents therefore 
very peculiar features. 
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Finally, wc can plot the magnetization 6 as a function 
of the energy e and polytropic inverse temperature rj by 
eliminating x between equations ([M)) . ([55)1 and ([55)) and 
between equations (jH]), (j45p and (pS)) . The corresponding 
curves are represented in Figures 0] and O According to 
the previous results, the inhomogeneous phase is unstable 
for 6 < ~ 0.37, metastable for < 6 < 5f ~ 0.44 and 
fully stable for b > bf. Note that the magnetization curves 
present a discontinuity at et and rjt confirming that the 
transition is first order. 



3.5.4 The control parameter e 

Let us summarize the previous discussion by expressing 
the dynamical stability of the waterbag distribution with 
respect to the Vlasov equation in terms of its energy e and 
magnetization b: 

(i) The spatially homogeneous waterbag distributions 
{b = 0) are fully stable for e > et — 0.352, metastable for 
€c = 1/3 <£<£(= 0.352 and unstable for e < Cc = 1/3. 

(ii) The spatially inhomogeneous waterbag distribu- 
tions are fully stable for emin = —2 < e < ft = 0.352 and 
b> bt= 0.44, metastable for et = 0.352 < e < = 0.379 
and = 0.37 < b < bt = 0.44. They do not exist for 
e > = 0.379. 

In principle, we cannot conclude that the spatially in- 
homogeneous waterbag distributions with ec ~ 1/3 < e < 
e* = 0.379 and 6 < 6* = 0.37 are Vlasov unstable since 
the criteria (|^5)) and P^)) provide just sufficient conditions 
of dynamical stability |47j . However, it seems natural that 
these solutions are unstable. 



3.5.5 The control parameter rjQ 




0.5 I 1.5 2 



Fig. 6. Energy e of the polytrope n = 1/2 as a function of 
the maximum value of the distribution function /i. Tliis corre- 
sponds to the ground state energy tground in the Lynden-Bell 
theory when the initial condition has only two levels and rio . 
The inhomogeneous branch tends towards e = —2 for ji — -\-oo. 
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Fig. 7. Zoom of Figure [6] near the bifurcation point. 



A polytrope of index n is characterized by its energy E (in 
dimensionless form e) and its inverse polytropic tempera- 
ture 1/K (in dimensionless form 77) . These are the control 
parameters that we have used in the previous sections. 
For the waterbag distribution, equivalent to a polytrope 
n = 1/2, it is better to use the maximum value of the 
distribution function 770 as a control parameter instead of 
the polytropic temperature K . Let us introduce the nor- 
malized maximum value of the distribution function 



27rfc 



1/2 



(67) 



According to equation it is related to the normalized 
inverse polytropic temperature (|25p by 



(68) 



We can now reformulate the preceding results in terms 
of the variables e and ^. The curve giving the energy e 
as a function of the maximum value of the distribution 
function ^ is plotted in Figures E] and [T] According to 
equations ([23 and ([55)) . the energy of the homogeneous 



phase is related to the maximum value of the distribution 
function by 



1 

6^' 



(69) 



The stability criterion (j60p expressed in terms of the max- 
imum value of the distribution function becomes 



^ S 



1 

72' 



(70) 



The homogeneous phase exists for ^ > and e > 0. It is 
fully stable for fi < nt — 0.688 (i.e. e > et), metastable 
for ^t < ^J' < ^J'c (i-C. ec < e < et) and unstable for 
/.t > i-ic (i.e. e < ec). The inhomogeneous phase exists 
for > fit 0.675 and —2 < e < e*. It bifurcates from 
the homogeneous phase at iic = l/\/2 and ec = 1/3. The 
upper part of the inhomogeneous branch is unstable. The 
lower part of the inhomogeneous branch is metastable for 
/I* < M Mt (i-C- e* <£<£*) and fully stable for fi > fit 
(i.e. emin < e < et)- The curve giving the magnetiza- 
tion 6 as a function of the maximum value of the distri- 
bution function fi is plotted in Figure [51 The inhomoge- 
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neous phase is unstable for b <b^, cri 0.37, metastable for 
6, < 6 < 6( ~ 0.44 and fully stable for h > bt. 
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Fig. 8. Magnetization 6 as a function of the maximum value 
of the distribution function fj,. 



Remark: recalling that the waterbag distribution (jl2[) 
corresponds to the minimum energy state (i.e. T = 0) 
in the two levels approximation of the Lynden-Bell the- 
ory [3HI I the curve of Figure El gives the ground state en- 
ergy Cground as a fuuctiou of the initial value ^ of the 
distribution function. As will be shown in the next sec- 
tion, the stable states (S) correspond to global minima 
of energy while the metastable states (M) correspond to 
local energy minima and the unstable states (U) to sad- 
dle points. In previous works |38l42l43j . only the min- 
imum energy eminifJ-) = 1/(6/^^) of the homogeneous 
phase, and the minimum energy eM/Af(A') of the inhomo- 
geneous phase in the particular case where the initial con- 
dition is a rectangular waterbag distribution, had been 
determined. Here, we have determined the absolute mini- 
mum energy egro«nd(M) of the inhomogeneous phase. For 
^ +(X) (non degenerate limit), we recover the classi- 
cal result CgroundilJ-) ^ ^ 2. Therefore, the curve of Fig- 
ure [S] completes the phase diagram of f43| by giving the 
minimum accessible energy (ground state). This point is 
specifically discussed in Appendix A of [44] . 

3.6 Thermodynamical stability of the Fermi 
distribution 

In the previous sections, we have studied the dynamical 
stability of the Fermi (or waterbag) distribution function 
and density profile ^7} with respect to the Vlasov 
and Euler equations. Wc shall now consider the thermo- 
dynamical stability of the Fermi distribution. 

3.6.1 The minimum energy state 

In the mean field approximation, the statistical equilib- 
rium state of a gas of fermions in interaction is obtained 
by maximizing the Fcrmi-Dirac entropy at fixed mass and 



energy in the microcanonical ensemble or by minimizing 
the Fermi-Dirac free energy at fixed mass in the canonical 
ensemble (see Section [O]) . Here, we are interested by the 
ground state E = E ground corresponding to T = 0. Wc 
thus have to determine the minimum energy state for a 
given value of mass while respecting the Pauli exclusion 
principle / < ?7o. 

To solve this minimization problem, wc can proceed 
in two steps. We first minimize the energy at fixed mass 
and density profile p{9) = J f{0, v) dv. Since the specifi- 
cation of the density profile determines the mass and the 
potential energy, this is equivalent to minimizing the ki- 
netic energy at fixed density profile. This is achieved by 
populating the lowest kinetic energy states with the maxi- 
mum of fermions allowed by the Pauli exclusion principle. 
This leads to the optimal distribution function 

/* = Vo, if \v\ < vpiP), 

=0, if \v\>VF{e), (71) 

where the Fermi velocity wf(^) is determined by the den- 
sity profile according to vp{9) = p{9)/{2r]o). We can now 
express the energy ([2]) as a functional of p by writing 
E[p] = E[f^]. This yields 

E^^-^ J p^de + ^J p<Pde, (72) 

where K is given by equation (jl9p . Finally, the ground 
state corresponds to the optimal distribution function 
with the optimal density profile that minimizes E[p\ at 
fixed mass M. 

Remark 1: we can proceed differently by considering 
the limit T — >■ 0+ of the Fermi-Dirac free energy. For T > 
0, the Fermi-Dirac free energy is given hy F = E — TS 
where E is given by equation ^ and S by equation (fTU|) . 
Now, for a general entropic functional of the form S = 
— J C{f) dOdv where C is convex, it is shown in |9l50l51j 
that the minimization of F[f] at fixed mass is equivalent 
to the minimization of F[p] at fixed mass, where F[p\ is 
the free energy (|5ip . In the case of fermions, p — p{p) 
is given by the Fermi-Dirac equation of state (see, e.g., 
[SS]). For r = 0, it reduces to the polytropic equation of 
state P?)l and the free energy F[p\ reduces to the energy 
functional ((72|l . 

Remark 2: for a polytropc of index n = 1/2, the free 
energy ([51]) is equivalent to the energy functional ([7^ . 
We can therefore directly use the results of Section 13.51 
to investigate the thermodynamical stability problem. In 
particular, the stability of the homogeneous phase is de- 
termined by the criteria given in Section 13.5.21 since they 
are precisely obtained by minimizing F[p] at fixed mass 
|9l47l54j . The stability of the inhomogeneous states can 
be determined by using the results of Section [3.5.31 or by 
proceeding as in the next section. 

Remark 3: a maximum of entropy at fixed mass and 
energy, a minimum of free energy at fixed mass and a min- 
imum of energy at fixed mass with the constraint f < rjQ 
are guaranteed to be dynamically stable with respect to 
the Euler and Vlasov equations. Therefore, thermody- 
namical stability implies dynamical stability. However, the 
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converse is not always true since the criteria (|3H1) and (|49|) 
provide just sufficient conditions of dynamical stability 

m- 



3.6.2 The series of equilibria 

The critical points of energy E, given by equation ([7^ . at 
fixed mass M are determined by the variational principle 
dE — a6M = where a is a Lagrange multiplier associ- 
ated with the conservation of mass. This returns the den- 
sity profile (P7)) with a = ep. In the following, we shall 
consider the case > since it corresponds to the region 
where the phase transition occurs. To determine whether a 
critical point of energy is a minimum or a saddle point, we 
can use the Poincare theorem. To that purpose, we have to 
plot a as a function of M (for fixed h) since a ~ dE/dM. 
Since a is an increasing function of A, this is equivalent 
to plotting A as a function of M (for fixed h). For the 
inhomogeneous phase, using the results of Section [XH we 
obtain the relations 



M 



A 



12k 



1 



12k 



(73) 



(74) 



which determine A{A4) in parametric form (with pa- 
rameter a;). For the homogeneous phase, we simply have 
A = Ai/2. The curve A{Ai) is represented in Figure HI 
We know from the study of Section 13.5.21 that the homo- 
geneous phase is stable for > Mc = 1/3 and unstable 
for < A^c- Then, using the Poincare theorem, we de- 
duce that the inhomogeneous phase is unstable close to 
the bifurcation point but that it becomes, and remains, 
stable after the turning point of mass at ~ 0.365. 



0.25 



0.2 - 




0.15 



Fig. 
M. 



9. Tlie chemical potential ^ as a function of the mass 



To study the stability of the Fermi distribution, we 
can also plot the energy £' as a function of the mass M 
and compare the energy of the different solutions that have 



0.05- 



0.045 - 



0.04- 



0.035 - 



0.32 0.33 0.34 0.35 0.36 0.37 
M 

Fig. 10. The energy £ as a function of the mass M. 



the same mass (for fixed h). For the inhomogeneous phase, 
using the results of Section 13.41 we obtain the relations 



£ = 



Eir^h* _ e(x) 



18P nixy 



M = 



12k 



-nix)'' 



(75) 



(76) 



which determines £{A4) in parametric form (with pa- 
rameter x). For the homogeneous phase, using equation 
we simply have £ = M^. The curve £{M) is rep- 
resented in Figure [101 This curve determines the transi- 
tion mass Ait — 0.352 separating stable and metastable 
states (here, fully stable states are global minima of en- 
ergy and metastable states are local minima of energy). 
The curves of Figures [HI and [TU] display a first order phase 
transition marked by the discontinuity of a{M) = E'{M) 
aX M = Mt and the occurrence of metastable states. On 
the other hand, can be interpreted as a spinodal point 
marking the end of the inhomogeneous metastable phase. 

In conclusion: the homogeneous phase exists for M > 
0. It is fully stable for M > Mt, metastable for Mc < 
M < Mt and unstable for M < Mc- The inhomoge- 
neous phase exists for < Al < A^,. It bifurcates from 
the homogeneous phase at M = Mc- The lower part of 
the inhomogeneous branch is unstable. The upper part of 
the inhomogeneous branch is fully stable ioi M < Mt 
and metastable for Mt < M < M*- Since M = 1/ry, 
we recover the same stability results as in the dynami- 
cal approach of Section 13.51 but from a different point of 
view. This equivalence was expected in view of Remark 2 
of Section 13.6.11 



3.6.3 The control parameter h 

In the previous section, we have fixed Ti and taken the 
mass M as a control parameter. This is the right way to 
solve the thcrmodynamical problem corresponding to the 
minimization of energy at fixed mass. Now that we have 
established the conditions of stability, in order to present 
the final results, it is more relevant to fix M and take ?i as a 
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control parameter. We therefore introduce the normalized 
Planck constant for fermions 



1 ' I ' I ' I ' I ' I ' r 



\ 2k 



X = ^ 



1/2 



(77) 



According to equations (|19p and it is related to the 
parameter ij by 



X 



1/2 



(78) 



Since M is fixed, the curves giving the energy and the mag- 
netization as a function of the normalized Planck constant 
correspond to e(x) and b{x)- They are plotted in Figures 
[TTl and [T21 We can now reformulate the preceding re- 
sults in terms of the normalized Planck constant x- Before 
that, we note that the normalized Planck constant is re- 
lated to the normalized maximum value of the distribution 
function by 



1 

X = -■ 



(79) 



Therefore, the results will coincide with those obtained in 
Section [333] up to a slight reinterpretation of the param- 
eters. However, we describe the results in detail in order 
to facilitate the comparison with the case of bosons in 
Section H31 



1.5 - 




0.5 - 



3.5 



Fig. 11. Ground state energy e as a function of the Planck 
constant X- The inhomogeneous branch (stable) tends towards 
e = — 2 for X — >■ 0. 



The curve giving the ground state energy e as a func- 
tion of the Planck constant x (in dimensionless units) is 
plotted in Figures [TTl and [T^ According to equations ([5^ 
and ([79]) . the ground state energy of the homogeneous 
phase is related to the Planck constant by 



X2 

6 



(80) 



On the other hand, the stability criterion ([5(7|) can be ex- 
pressed in terms of the Planck constant as 



(Z.,E.) 




Fig. 12. Zoom of Figure [TT] near the bifurcation point. 



0.8 



0.6 



0.4 



0.2 



1 









1 




~~— — — 





















M 








^ (X.,b.) - 














U / 






U 


X, 1 1 M 


Z, S 
1 


2 


1.3 


1.4 


1.5 



X 



(81) 



Fig. 13. Magnetization & as a function of the Planck constant 
X- It displays a first order phase transition between the classical 
regime (inhomogeneous state & 7^ 0) and the quantum regime 
(homogeneous state b — 0). 



The homogeneous phase exists for x > and e > 0. It is 
unstable for x < Xc (i-S- e < €c), metastable for Xc < X ^ 
Xt — 1-45 (i.e. Cc < e < Ci ~ 0.352) and fully stable for 
X > Xt (i-C. e > 64). The inhomogeneous phase exists for 
X < X* — 1-48 and — 2 < e < e* ~ 0.379. It bifurcates from 
the homogeneous phase at Xc = and ec — 1/3. The 
upper part of the inhomogeneous branch is unstable. The 
lower part of the inhomogeneous branch is metastable for 
Xt < X ^ X* (i-C- and fully stable for < x < 

Xt (i.e. —2 < e < et). The curve giving the magnetization 
b versus the Planck constant x is plotted in Figure [T31 
The inhomogeneous phase is unstable for 6 < ~ 0.37, 
metastable for 6* < & < &t ~ 0.44 and fully stable for 
6 > It displays a first order phase transition at x = Xt 
marked by the discontinuity of the magnetization. 

Summarizing, we have determined the ground state 
(T = 0) of the Fermi gas with cosine interaction and inves- 
tigated its thermodynamical stability. In the classical limit 
X = 0, the ground state energy is e = — 2 and the magneti- 
zation is & = 1. This corresponds to a classical inhomoge- 
neous gas at T = whose density profile is a Dirac peak 
p = M5{9). The homogeneous phase with energy e = 
and magnetization 6 = is unstable since T = < Tc. 
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This returns the classical results [3]. However, we note 
that when quantum mechanics is taken into account, there 
exists a critical value of the normalized Planck constant 
above which the homogeneous phase becomes stable (more 
precisely, it is metastable for Xc < X < Xt and fully stable 
for X > Xt). In parallel, the inhomogeneous phase becomes 
metastable for Xt < X < X* and disappears for x > X*- 
Therefore, in the quantum regime, the homogeneous phase 
is stabilized against clustering by the Pauli exclusion prin- 
ciple. This is similar to the stabilization of fermion stars 
against gravitational collapse in astrophysics due to quan- 
tum mechanics. In the HMF model, this stabilization takes 
place through a first order phase transition. 

Remark: the various approaches developed in Sections 
l3.5l and l3.6l show the importance of the parameter 77 in the 
stability analysis. This parameter can have several inter- 
pretations: inverse polytropic temperature maximum 
value of the initial condition inverse mass ([75| and 
inverse Planck constant ([75)1 . The change of stability cor- 
responds to the turning point of rjix) that takes place at 
77* ~ 2.737. This is relatively similar to results obtained 
in astrophysics [55] . 



BECs |30l33l34j . The validity of the Schrodinger-Poisson 
system has been justified rigorously by Lieb & Yau [57] . 
In this paper, we shall study the bosonic HMF model at 
T — Q that is described by the mean field Schrodinger 
equation ([5^ with a cosine potential of interaction (jS]). 



4.2 The Madelung transformation 

Let use the Madelung jSS] transformation to rewrite the 
Schrodinger equation in the form of hydrodynamic equa- 
tions. We first set 



(86) 



where A(r, t) and ^(r, t) are real functions. We clearly 
have 



A 



2i 



(87) 



where 1}}* denotes the complex conjugate. Substituting 
equation ()86p in equation (j82p and separating real and 
imaginary parts, we obtain 



4 Bosons at zero temperature 

In this section, we consider the HMF model at T = in 
the case where the particles are bosons. We use a mean 
field approximation that becomes exact in a proper ther- 
modynamic limit N — > +CX) defined in Appendix El 



4.1 The mean field Schrodinger equation 

Let us consider a system of bosons interacting via a long- 
range binary potential w(|r — r'|). At zero temperature 
T = 0, these particles form a Bose-Einstein condensate 
(BEC) described by a wave function 7/;(r, t). For N +00, 
the wave function is solution of the mean field Schrodinger 
equation 

in-^{r,t) = -—Ai;{r,t)+m^{v,t)ij{r,t), (82) 



<?(r,0 = / pir',t)ui\r-v'\)dr', 



pir,t) = Nm\ij{r,t)f' 



(83) 



(84) 



dS_ 

dt 



1 
2m 

dt 



{VSf + m$ 



f^AA 
2m~~A 



= 0, 



0. 



(88) 



(89) 



The first equation has a form similar to the classical 
Hamilton- Jacobi equation with an additional potential 
term Q = ^ 5^ called the quantum potential. Follow- 
ing Madelung, we introduce the density and velocity fields 



p = NmA^ = Nm\iPf 



u = —VS. 

m 



(90) 



We note that the flow defined in this way is irrotational 
since V x u = 0. With these notations, equation ([5^1) 
becomes similar to the equation of continuity in hydrody- 
namics. On the other hand, equation ([55)1 can be inter- 
preted as a generalized Bcrnouilli equation for a potential 
flow. Taking the gradient of equation ((88)) and using the 
well-known identity (u • V)u = V(^) — ux (Vxu) which 

reduces to (u- V)u = V(^) for an irrotational flow, we 
obtain an equation similar to the pressureless Euler equa- 
tion with an additional quantum potential. In conclusion, 
the mean field Schrodinger equation is equivalent to the 
"hydrodynamic" equations 



mr,t)\Ur = l. 



(85) 



Equation (j85p is the normalization condition, equation 
([M)) gives the density profile of the BEC, equation (|83p 
determines the associated potential and equation ([5^ de- 
termines the wavefunction. For the gravitational poten- 
tial, these equations reduce to the Schrodinger-Poisson 
system which describes boson stars and self-gravitating 



where 

Q 



|+V.(pu)^0, 



^ + (u.V)u 



Am 



2m yfp 



-V<?- -VQ, 

TO 



A^ 
P 



nvpf 

2 p2 



(91) 



(92) 



(93) 
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is the quantum potential. We shall refer to equations 
((92|) and ((83)) as the quantum Euler equations. In the clas- 
sical limit ?i — ^ 0, the quantum potential disappears and 
we recover the ordinary Euler equations. 

Remark: the quantum potential first appeared in the 
work of Madelung [S^ and was rediscovered by Bohm [SH] 
(it is sometimes called the Bohm potential). We note the 
identity 

^-VQ = --d,P,,, (94) 



where is the quantum stress tensor 



p>. . — 



pdidj Inp, 



or 



(95) 



(96) 



This shows that the quantum potential is equivalent to an 
anisotropic pressure. 



4.3 The time independent Schrodinger equation 

If we consider a wavefunction of the form 

4>{r,t) = A{r)e-^^ , (97) 
we obtain the time independent Schrodinger equation 



2m 



(98) 



where ip{r) = ^(r) is real and p{r) = Nrmjj'^iY). Dividing 
equation ((55)1 by ■(/'(r), we get 



AJp 



(99) 



which can be written in terms of the quantum potential 
as 

m^ + Q = E. (100) 

This relation can also be derived from the Hamilton- 
Jacobi equation ((55)) by setting 5* = —Et. 

Using the Madelung transformation, equation pOOp 
represents the steady state of the quantum Euler equa- 
tions ((OT)) -((M | . Indeed, taking (9f = and u = in equa- 
tion ((HU, we get 



V<?+ -VQ = 0, 



(101) 



which can be interpreted as a condition of quantum hy- 
drostatic equilibrium. It describes the balance between the 
long-range interaction and the quantum pressure due to 
the Heisenberg uncertainty principle. Integrating this re- 
lation, we recover equation ()100p where the energy E ap- 
pears as a constant of integration. 



4.4 The total energy 

The energy functional associated with the mean field 
Schrodinger equation, or equivalently with the quantum 
Euler equations, is 

Etot ^Oc + Oq + W. (102) 

The first two terms correspond to the total kinetic energy 



2m 



(103) 



Using the Madelung transformation, it can be decomposed 
into the "classical" kinetic energy 



u 



0c = / P Y dr. 



and the "quantum" kinetic energy 

1 

TO 



0C 



pQ dr. 



(104) 



(105) 



Substituting the expression of the quantum potential ((931 
in equation ()105p . we equivalently have 



Q 



2to2 



y/pA^dr 



2to2 ./ 8to2 
The third term is the potential energy 

W ^ 



{Vpf 



p(l>dr. 



dr. (106) 



(107) 



The total energy can be expressed in terms of the wave- 
function as 



Etnt — 



^\^^\2 + ^Nmm^ 
Zm 2 



(108) 



and the mean field Schrodinger equation can be written 



ih- 



dt 



(109) 



It is easy to show that the energy functional ()102p and 
the total mass M = J pdr are conserved by the quantum 
Euler equations. This implies that a minimum of Etot at 
fixed mass M is a nonlinearly dynamically stable steady 
state of the quantum Euler equations [48]. Writing the 
first variations as SEtot ~ aSM = 0, where a is a Lagrange 
multiplier, we get u = and 



m(p + Q = ma. 



(110) 



Taking the gradient of this relation, we obtain equation 
((lOip which characterizes a steady state of the quantum 
Euler equations. Equation ()110p also coincides with the 
time independent Schrodinger equation (|100p provided 
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that the Lagrange multipUer a and the energy E are re- 
lated to each other by E = ma. Therefore, the energy 
E arising in the time-independent Schrodinger equation 
can be interpreted as a chemical potential. At equilibrium 
(u = 0, 0c = 0), the total energy reduces to 



Q 



W. 



(Ill) 



Multiplying equation (|100p by p and integrating over the 
entire domain, we obtain 



Therefore, 



2W + eQ^ NE. 



Etot ^NE-W. 



(112) 



(113) 



From this equation we note that Etot 7^ NE when W ^ 0. 
In the classical limit h ^ 0, 0q = 0, so that 



Etot = = l-NE. (classical) 



(114) 



4.5 The cosine interaction 



For the bosonic HMF model at T = 0, corresponding 
to a cosine interaction, the time independent Schrodinger 
equation (|98| takes the form 



-ij"{e) + ^e)i,{e) = E^{0), 



(115) 



where 



m 



27r Jo 



p{9') cos{0 - e')d9', (116) 



(117) 



2ir 



V'^ d9 = 1, 



Pi9) - Ni'^9), 



(123) 



(124) 



where b = 2TiB/kM is the magnetization. The boundary 
conditions are 



^'(0) = ^'(Tr) = 0. 



(125) 



This system of equations defines an eigenvalue problem 
for the wave function ■0(^') where the eigenvalue E is the 
energy. In the following, we shall be interested by the fun- 
damental eigenmode corresponding to the smallest value 
of E. For this mode, the wave function "0(0) has no node 
so that the density profile decreases monotonically with 
the distance. 




-0.6 



Fig. 14. Eigenvalue A as a function of k^. The dashed line 
starting at = corresponds to the harmonic approximation 
of Appendix [P] 



Jo 



ip^{9)d9 



(118) 



We can suppose without loss of generality that ip{9) and 
p(9) are symmetric with respect to 9 = 0. In that case, 
'P{9) is given by 



m 



-B cosl 



with 



2lT 



B = — / p{9) COS 9 d9. 



(119) 



(120) 



Therefore, the time independent Schrodinger equation 
with a cosine potential can be rewritten 



—V'" - —bcos9tP Eij, 
2 27r 



2ir 



COS 9 d9, 



(121) 
(122) 




Fig. 15. Magnetization 6 as a function of k^. 



We note that ^p{9) = l/v27r is a particular solution of 
these equations. It corresponds to a spatially homogeneous 
distribution p{9) = M/{2tt) with 6 = 0. Its stability 
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Fig. 16. Normalized Planck constant x as a function of k . The 
normalized Planck constant tends to the critical value Xc ~ \/2 
when K — s> oo. 




the fundamental eigenvalue A = A(k) corresponding to the 
wavefunction that has no node. Then, we can obtain the 
magnetization from the equation 



/•27r 

b= iP^cosOde = bin). 
Jo 



Finally, we get 



\kM J 



2ttE 



1/2 



kM 



A(k)6(k). 



(130) 



(131) 



(132) 



These equations determines the normalized energy £ = 
£{x) and the magnetization b = b{x) as a function of the 
normalized Planck constant x in parametric form (with 
parameter k). In fact, we can simplify the problem even 
further. Indeed, if we define (f>{0) = ip{d)/ip{0), we obtain 



(133) 



-(/)" — cos 9(j) — X(j), 



0(0) = 1, 0'(O) = 0'(7r) = 0. 



(134) 



Fig. 17. Normalized energy £ as a function of 



This is a simple shooting problem. For given k, solving 
equation p33p with the initial condition 4>{0) ~ 1 and 
0'(O) = 0, we have to find the eigenvalue A(k) so as to 
satisfy the boundary condition (/''(tt) = 0. Then, ^/'(O) is 
given by the normalization condition (|129|) leading to 



is investigated in Appendix [Bl Introducing the normalized 
Planck constant for bosons 



1/2 



x = h 



(126) 



\kM J ' 

we find that the homogeneous phase is stable if 

X>Xc = V2, (127) 

and unstable otherwise. We stress that the normalized 
Planck constant is different for fermions and bosons (see 
Appendix 

Inhomogencous solutions with b must be obtained 
numerically. Dividing equation (|12ip by fcM6/(27r) and 
setting 2nh'^/{kMb) and A = 2nE/{kMb), we obtain 



— V'" — cosOip = Xip, 



ijj^ de = 1. 



(128) 



(129) 



For given k, this is just the ordinary Schrodinger equation 
for an anharmonic oscillator. We can therefore determine 



(135) 



This completely determines the wave function ijj{9) = 
ip{Q)(l){d) and the density profile p{e)/N = i^iOf . Finally, 
the magnetization is given by 



(136) 



while the normalized Planck constant and the normalized 
energy are given by equations (|13ip and (|132p . 

We recall that the energy E appearing in the time inde- 
pendent Schrodinger equation ()115p is generally different 
from the total energy per particle Etot /N (see Sectional]). 
The total energy is given by equation (jll3|) . For the HMF 
model, the potential energy can be expressed in terms of 
the magnetization by equation ([7]). Therefore, the total 
energy is 

Etot = NE+—. (137) 
Introducing the normalized energy 



Etot 



kM^ ' 



(138) 
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Fig. 18. Magnetization as a function of the normalized 
Planck constant x- The dashed line starting at x = cor- 
responds to the harmonic approximation of Appendix |D] The 
dashed line starting at x = Xc = \/2 corresponds to the an- 
alytical expression of the magnetization close to the critical 
point (see Appendix [C]). The magnetized solution exists only 
for X < Xc and it is stable. The homogeneous phase (6 = 0) is 
unstable for X < Xc and stable for X > Xc (see Appendix [B)) . 
This curve displays a second order phase transition between the 
classical regime (inhomogeneous state fe 7^ 0) and the quantum 
regime (homogeneous state & = 0). It can be compared with 
the magnetization curve of fermions reported in Figure [THl 




Fig. 19. Normalized energy £ as a function of the normalized 
Planck constant x- 



we get 

ttot = 4f + 2h\ (139) 

In the homogeneous phase where 6 = f = 0, we have 
Ctot = for any value of x- This differs from the case of 
fermions where the energy of the homogeneous phase is 
given by equation (1501) . 

The curves A(k), 6(k), x(k) and E{k) are plotted in 
Figures [H HH [11 and [HI The curves 6(x), f (x) and 
e*ot (x) giving the magnetization and the normalized ener- 
gies as a function of the normalized Planck constant are 
plotted in Figures [TSl [10] and [501 Finally, some density 
profiles are shown in Figure [2TJ The homogeneous phase 
exists for any x- It is stable for x > Xc = and unstable 
for X < Xc (sec Appendix |B|) . The inhomogeneous phase 




jk ^ \ ^ \ ^ \ ^ I 

0.5 I 1.5 2 

X 

Fig. 20. Total normalized energy etot as a function of the 
normalized Planck constant x- It can be compared with the 
energy curve of fermions reported in Figure 1111 




Fig. 21. Density profiles for x = 0.0376, 0.1, 0.3, 0.85, 1.33, ^/2 
(from top to bottom). The dashed lines correspond to the 
Gaussian approximation of Appendix [D] These profiles can 
be compared with the density profiles of fermions reported in 
Figure [T] 



exists only for x < Xc = where it bifurcates from the 
homogeneous phase. The behaviors of the magnetization 
and of the energy close to the bifurcation point are ob- 
tained analytically in Appendix [C] and compared with the 
numerical results in Figures [15] and [TO] Analytical results 
are also obtained in the semi-classical limit x ~^ by mak- 
ing a harmonic approximation (see Appendix |D]) . These 
analytical results are compared with the numerical ones 
in Figures [T51I?I1 According to the Poincare theorem, we 
deduce that the inhomogeneous phase becomes stable at 
the point X = Xc — \/2 at which the homogeneous branch 
becomes unstable. Since the inhomogeneous branch is not 
multivalued, there is no metastable state contrary to the 
case of fermions. The system displays a second order phase 
transition at x X* marked by the discontinuity of the 
derivative of the magnetization. 

Summarizing, we have determined the ground state 
(T = 0) of the Bosc gas with cosine interaction, forming 
a BEC, and investigated its thermodynamical stability. 
In the classical limit x = 0, the ground state energy is 
Etot = — 2 and the magnetization is 6 = 1. This corre- 
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sponds to a classical inhomogcneous gas at T = whose 
density profile is a Dirac peak p = M5{9). The homoge- 
neous phase with energy ttot = and magnetization 6 = 
is unstable since T = < Tj,. This returns the classical 
results [5]. However, when quantum mechanics is taken 
into account, there exists a critical value of the normal- 
ized Planck constant Xc = above which the homo- 
geneous phase becomes stable while the inhomogcneous 
phase disappears. Therefore, in the quantum regime, the 
homogeneous phase is stabilized against clustering by the 
Heisenberg uncertainty principle. This is similar to the 
stabilization of boson stars against gravitational collapse 
in astrophysics due to quantum mechanics. In the HMF 
model, this stabilization takes place through a second or- 
der phase transition. 



p{p) is determined by the potential h{p) through the rela- 
tion 

h\p) - (146) 
P 

yielding p{p) = ph{p) — H{p) where H is a primitive of 

h. Using this relation, the Euler equation ()145|) can be 

rewritten as 

o -1 

^ -f (u • V)u = -V/i - V<? VQ, (147) 

at m 

which shows that the potential h appearing in the GP 
equation can be interpreted as an enthalpy in the hydro- 
dynamic equations. Equations (|144|) . (|145p and (|14ip will 
be called the quantum barotropic Euler equations. They 
are equivalent to the mean field GP equation. 



5 The mean field Gross-Pitaevskii equation 

We shall now consider a generalization of the equations of 
Section m by introducing a nonlinearity in the Schrodinger 
equation in addition to the long-range interaction. This is 
the so-called mean field Gross-Pitaevskii (GP) equation. 
As we shall see, this equation can have application for 
fermions and bosons. For the sake of generality, we con- 
sider an arbitrary potential of interaction in d dimensions. 

5.1 The Madelung transformation 

The mean field Gross-Pitaevskii equation can be written 
in^ = -^M' + m{^ + h{p))2p, (140) 



<P{r,t) = / p{r',t)u{\r^r'\)dr', (141) 



p{r,t) = Nm\^{r,t)\^, (142) 

y|7/.(r,t)|2rfr = l, (143) 

where h(p) is a potential depending on the density and the 
other quantities have been defined previously. Adapting 
the Madelung transformation of Section HT^ to the present 
situation, we find that the Gross-Pitaevskii equation is 
equivalent to the hydrodynamic equations 

|^-fV.(pu)=0, (144) 

O I -1 

^ + (u-V)u = — Vp-V^ VQ, (145) 

ot p m 

where Q is the quantum potential (|93| and p(r, t) is a pres- 
sure that is a function of the density: p{r,t) = p[p{r,t)]. 
In this sense, the fluid is barotropic. The equation of state 



5.2 Interpretations of the potential 

We can give two independent interpretations of the po- 
tential h(p). 

(i) Let us consider a BEG described by the mean field 
Schrodinger equation (|5^ - ((55|) and let us assume that the 
potential of interaction can be written a.s u = ulb. + usr 
where ulr refers to the long-range interaction and usr 
to the short-range interaction. We assume furthermore 
that the short-range interaction corresponds to binary 
collisions that can be modeled by the effective potential 
usr{^ — r') = gS{r — r'), where the coupling constant (or 
pseudopotcntial) g is related to the scattering length a by 
g = iTTati^/mP (in d = 3) [5D]. When this form of poten- 
tial is substituted in equation (j83p . we obtain the usual 
mean field Gross-Pitaevskii equation 

^n^ = -^A^ + m'Pij + N^^^\ij\''ij. (148) 
ot 2m m 

with a potential h = gp = gNm\ip\'^. The associated equa- 
tion of state is p = ^gp^ corresponding to a polytrope of 
index n = 1 and polytropic constant K = g/2. 

(ii) The mean field Gross-Pitaevskii equation (|140p - 
(|143p also describes a gas of fermionf[£| when one takes into 
account the quantum potential Q which is a manifestation 
of the Heisenberg uncertainty principle |61j . In the case of 
fermions, we must also take into account the quantum 
pressure arising from the Pauli exclusion principle. It can 
be calculated from the Fermi-Dirac distribution function 
at r = (see Section In d dimensions, it is given 
hyp = Kp^+y'^ where K = ^^{^)Vd^0L ggj. This 
is the equation of state of a polytrope of index n = d/2 
and polytropic constant K. This pressure term is the one 
that appears in the hydrodynamic equation (|145p . Using 

^ Such a description assumes that the fermions have the same 
probability distribution [ST]. More fundamentally, the fermions 
should be described by a mixture of A'^ pure states, each with a 
wavefunction obeying the mean field Schrodinger equation 
without nonlinearity (except the one due to the interaction) 
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equation (jl46p . it corresponds to an effective potential of 
the form h{p) = {d/2+l)Kp^^'^. The corresponding Gross- 
Pitaevskii equation can be written 




where Kd = 2tt^{^-^)'^^'^ is a constant. 

Remark: we note that the potential h{p) cx p^^'^ as- 
sociated to fermions becomes equivalent to the potential 
h{p) oc p associated to self-coupled bosons when d = 2. 
In fact, the dimension d = 2 is a critical dimension 
[55] . When we consider a gas of repulsive (impenetrable) 
bosons, the potential h(p) = gp arising in the GP equa- 
tion ceases to be valid for d < 2 (in c? = 2 it remains 
marginally valid with logarithmic corrections). In partic- 
ular, in d = 1, it is replaced by h{p) ~ n^h p^ /2m'^ exactly 
like for spinless (s = 0) fermions. This is a manifestation 
of the boson- fermion duality in one dimension [64] . 

5.3 The time independent Gross-Pitaevskii equation 

If we consider a wavefunction of the form 

ip{r,t) = A{r)e-^^, (150) 
we obtain the time independent GP equation 

- —A^P{r) + m(<P(r) + h{p))iP{r) = EiP{r), (151) 
2m 

where ip{r) = A{r) is real and p{r) = Nrail)'^{Y). Dividing 
equation (|15ip by i/'(r), we get 

m<P + mh{p) - = E. (152) 
2m yfp 

This equation can also be obtained from the quantum 
barotropic Euler equations p44p . (jl45p and (|14ip since 
they are equivalent to the GP equation. The steady states 
of the Euler equation (|145p obtained by taking 9f = and 
u = satisfy 

VP + PV<.-|^V(^)=0. (153) 

This is similar to the condition of hydrostatic equilib- 
rium with an additional quantum potential (as we have 
previously indicated, it can be written in the form of an 
anisotropic pressure) . This equation is equivalent to equa- 
tion (|152p . Indeed, integrating equation (|153p and using 
equation (jl46p , we obtain equation (|152p where the energy 
E appears as a constant of integration. 

The equations of the problem can be simplified in two 
important limits: 

(i) The non interacting limit: the case of a BEC with- 
out short-range interactions corresponds to h = p = 
leading to 



This is the situation considered in Sectional In that case, 
the equilibrium state results from a balance between the 
long-range potential and the quantum pressure arising 
from the Heisenberg principle. 

(ii) The Thomas-Fermi limit: if we neglect the quan- 
tum potential in equation (|153p . we obtain 

Vp + pV<? = 0. (155) 

This corresponds to the standard condition of hydrostatic 
equilibrium describing a balance between the long-range 
potential and the pressure due (i) to the short-range in- 
teraction in a BEC (ii) to the quantum pressure resulting 
from the Pauli exclusion principle for fermions. 

Remark: considering the generalized condition of hy- 
drostatic equilibrium (jl53p taking into account the quan- 
tum potential, the Thomas-Fermi approximation is valid 
when the first term dominates over the third one. In the 
case of fermions, the pressure is due to the Pauli exclusion 
principle. Using the results of Section [5.2l and dimensional 
analysis, we easily see that the quantum potential (Heisen- 
berg uncertainty principle) can be neglected in front of the 
Fermi pressure (Pauli exclusion principle) when A'^ ^ 1. 
Therefore, the Thomas- Fermi approximation is exact for 
fermions in the thermodynamic limit N — > +oo. This is 
precisely the situation considered in Section [31 

5.4 The total energy 

The energy functional associated with the quantum 
barotropic Euler equations (jl44p . (jl45p and p4ip is 

Etot ^0c + Oq + U + W, (156) 

where 

U = j pj'P^dp,dv 

- / [ph{p)-p{p)] dr^J H{p)dv. (157) 

is the internal energy and the other functionals have al- 
ready been defined in Section 14.41 For a polytropic equa- 
tion of state p = Kp'^ , corresponding to a potential 
h{p) = ■^ZiP^^^i the internal energy takes the form 

U = I p'' dv = I pdr. (158) 
7 - 1 J 7-17 

It is easy to check that the total energy p56p and the 
mass AI ~ J pdv are conserved by the quantum barotropic 
Euler equations (hence by the GP equation). This implies 
that a minimum of Etot at fixed mass M is a nonlinearly 
dynamically stable steady state of the mean field quantum 
barotropic Euler equations [IS]- Writing the first varia- 
tions as SEfot — ol5M = 0, where a is a Lagrange multi- 
plier, we get u = and 

m<^ + mh{p) — = ma. (159) 

2m y/p 
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Taking the gradient of this expression and using equa- 
tion ()146|) we obtain the condition of hydrostatic equi- 
hbrium (|153p determining a steady state of the quantum 
barotropic Euler equations. This equation is also equiva- 
lent to the time independent GP equation (jl52p . The La- 
grange multiplier a is related to the energy E appearing 
in the time independent GP equation (jl52|) by a = E/m. 
This shows that the Lagrange multiplier a is equal to the 
energy E by unit of mass. On the other hand, consider- 
ing the second order variations of energy, we find that the 
distribution is dynamically stable iff 



quantum potential can be neglected and we recover the 
classical barotropic Euler equations considered in |9I53| . 

Linearizing equations (|165p - (|167p around a spatially 
homogeneous steady state (p = est and u = 0), we obtain 



9i5u h 



(168) 



(169) 



S^E,. 



dr 



h'{p){Spfdv 



„3/2 



Sp6<P dr 



S^{r,t) = j u{\v-v'\)5p{r',t)dv', (170) 



8m2 



{5pY dr > 0, where we have introduced the velocity of sound 



(160) 



for all perturbations that conserve mass: J Spdr = 0. 

At equilibrium (u = 0,0c — 0), the total energy re- 
duces to 

Etot = eQ + U + W. (161) 

Let us consider the case of a polytropic equation of state 
for which the internal energy is given by equation (jl58p . 
Multiplying equation (|152|) by p and integrating over the 
entire domain, we obtain 



-/U + 2W + 0Q= NE. 



(162) 



In the TF limit {Oq = 0), equations and be- 
come 

Etot ^U + W, (163) 

-/U + 2W = NE, (164) 

and in the noninteracting limit, we recover equations 
(fTTT]l-(ITT2]). 



5.5 Stability of the homogeneous phase with respect 
to the quantum Euler equations 

Let us study the linear dynamical stability of a spatially 
homogeneous distribution with respect to the quantum 
barotropic Euler equations 



|+V.(pu)=0, 



(165) 



^ + (u.V)u 



-Vp - p\/<P 



ph 



2™2 \ ^p ) 



(166) 



<P{r,t)= J u{\r~r'\)p{r',t)dr'. (167) 

These equations are equivalent to the mean field GP equa- 
tions (|Tig)) - (fTi5)) . In the classical or TF hmit h ^ 0, the 



cl=p'{p)^ph'{p). 
Equations (fT68)) and (fT69)) can be combined to give 



(171) 



d^p 



= clASp + pAS^ - 



Aim? 



A^5p. 



(172) 



Decomposing the perturbations in normal modes 5p 



uj'^Sp = -clk^Sp - pk^Sfp 



5$ = i2TTfu{k)dp. 
This leads to the dispersion relation 

a;2 = c^fc2 + {2Trf u{k)pk^ - 



Am? 



k'^6p, (173) 



(174) 



Am? 



(175) 



In the classical or TF limit ?i — >■ with Cg fixed, we recover 
the classical dispersion relation associated with the lin- 
earized Euler equation |9I53) . We note that the quantum 
pressure arising from the Heisenberg uncertainty principle 
(last term in equation (|175p ) becomes important at very 
small scales, i.e. k — >■ -|-cxd. According to equation (|175|) . 
the pulsation u is either real {u? > 0) or purely imaginary 
(w^ < 0). As a result, the system is linearly stable with 
respect to a perturbation with wavcnumbcr k if 



c? + (27r)'^H(fc)p+0>O, 



(176) 



and linearly unstable otherwise (when > 0, the pertur- 
bation oscillates with a pulsation uj and when o;^ < 0, the 
perturbation grows exponentially rapidly with a growth 
rate 7 = V— w^). 

For repulsive potentials, satisfying u{k) > 0, the ho- 
mogeneous distribution is always stable. For Coulombian 
plasmas in d = 3 dimensions, using {2TT)^u{k) = , the 
dispersion relation can be written 



,2^2 



ctk 



Am^ 



(177) 
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where uj^ = Airpe^ /m? is the plasma pulsation. This cor- 
responds to the Bogoliubov energy spectrum of the exci- 
tation of the boson ground state |65j . For large wavenum- 
bers (small wavelengths), the quasi-particle energy tends 
to the kinetic energy of an individual gas particle and 
uj ~ hk^/{2m). For the repulsive HMF model, using 



-Sn,±i, the dispersion relation can be written 



2 2 2 I ^-^'^ 2 jr 

UJ = c^n + — — n On,±i 
47r 



2^4 



The modes n 7^ ±1 oscillate with a pulsation uP' 



(178) 



h n /4. The pulsation of the mode n = ±1 is given by 

is the proper pulsation 



^ where uj^ 



kM 
i-rr 



For attractive potentials, satisfying u{k) < 0, the ho- 
mogeneous distribution is linearly stable if 



cl > {cl)crit = max 



{27rf\u{k)\p- 



(179) 



and linearly unstable otherwise. In that case, the unstable 
wavelengths are determined by the converse of inequality 
(jl76p . For the gravitational interaction in d = 3 dimen- 
sions, using {2Tr)^u{k) 
can be written 



the dispersion relation 



2 1,4: 



n^k 

Am? 



(180) 



The system is always unstable {{cDcrit = 00) for 
wavenumbers k < k^ where fc* is the generalized Jeans 
wavenumber [34l: 



and unstable otherwise. We note that quantum mechan- 
ics favors the stability of the homogeneous phase. If we 
consider a BEC without short-range interaction {h = p = 
Cs — 0), the stability criterion becomes 



kM^ tt' 



(186) 



and we recover the result of Appendix [Bl On the other 
hand, in the classical or TF approximation 7i — >■ with Cs 
fixed, the stability criterion becomes 



kM 
1^' 



(187) 



and we recover the result of [5]. 

Let us specifically discuss the case of fermions at T = 
described by the Fermi distribution ([T^ leading to the 
polytropic equation of state ((TT)) . The velocity of sound is 



Mh 

4 ' 



(188) 



and the pulsation equation for the modes n = ±1 can be 
written 



OJ = Vp - 



kM 
An 



A 



(189) 



The first term on the right hand side corresponds to the 
Fermi pressure (Pauli exclusion principle), the second to 
the cosine attraction and the third term to the quantum 
pressure (Heisenberg uncertainty principle) . In the TF ap- 
proximation, the pulsation equation reduces to 



2m2 



AnGhp 



(181) 



In the classical or TF limit ?i — > with fixed, we recover 
the classical Jeans wavenumber 



(182) 



For a non interacting BEC (cg =0), we obtain the quan- 
tum Jeans wavenumber 



kn = 



^IGnGm^p 



1/4 



(183) 



For the attractive HMF model, using ti„ ~ —^Sn,±i, the 
dispersion relation is 



2 22 
UJ = cji 



kM 2 c ^^n^ 

— n dn.±l -I — . 

47r 4 



(184) 



The modes rt 7^ ±1 oscillate with a pulsation uJ^ = c?^r? + 
h n^jA. The complex pulsation of the mode n = ±1 is 
uj^ ~ — + The homogeneous phase is stable for 



kM 
Att 



A ' 



(185) 



2 2 kM 
and the stability criterion (|187p can be written 



Mh^ A 

k TT 



(190) 



(191) 



This returns the result of Section [3] If we now take the 
quantum potential into account and consider the complete 
stability criterion (|185p . we obtain 



Mh' 



> 



(192) 



Comparing equation (jl92p with equation (jl9ip . we see 
that the term coming from the quantum potential has a 
contribution of order 1/7V. Since the mean field approx- 
imation on which our approach is based assumes that 
N — >■ -|-oo, it is not consistent to keep terms of order 
1/A^ in the calculations (other terms of order 1/A^ may 
appear as a deviation from the mean field limit). There- 
fore, the quantum potential can be neglected in the limit 
— > -l-oo. This is a general result. In the case of fermions, 
the quantum pressure arising from the Heisenberg prin- 
ciple is always negligible with respect to the quantum 
pressure arising from the Pauli principle for N +00 
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(see Section l5.3p . Therefore, the Thomas- Fermi approxi- 
mation of Seetion[3]is exact in the proper thermodynamic 
hmit defined in Appendix [Xj Alternatively, in the case of 
bosons without short-range interaction, the quantum po- 
tential stabilizes the system at small scales and must be 
taken into account in the thermodynamic limit defined in 
Appendix [21 This is the reason why fermions and bosons 
have different thermodynamic limits. 

5.6 Stability of the homogeneous phase with respect 
to the quantum Vlasov equation 

To complete our analysis, we study the linear dynami- 
cal stability of a spatially homogeneous distribution func- 
tion / = /(v) with respect to the Wigner equation. The 
Wigner approach is a reformulation of quantum mechanics 
in the classical phase space language. We recall, however, 
that Wigner functions can take negative values so that 
they are not true probability distributions. In the case of 
fermions, the Wigner equation, which can be viewed as a 
quantum Vlasov equation, is more fundamental than the 
quantum Euler equation considered in Section 15.51 which 
is based on approximations (see Section [5T2|) . For simplic- 
ity, we restrict ourselves to the one dimensional case. The 
Wigner equation reads 



dt dx 2nfi 



<P X 



Xh \ ( Xh 



<Pix,t) ^ u{\x-x'\)p{x',t)dx'. 



(193) 



(194) 



We can check that, in the classical limit h 0, equation 
(|193p returns the usual Vlasov equation. The dispersion 
relation corresponding to the linearized Wigner equation 
is (Ml: 



e{k,uj) EE 1 - 27ru(fc) 



f{v)dv 



4m^ 



0. (195) 



For ?i — > 0, we recover the classical dispersion relation 
associated with the linearized Vlasov equation j9l53j . The 
point of marginal stability (w = 0) is determined by the 
condition 



1 - 2TTu{k) 



f{v)dv 

2 _ h^k^ 



0. 



(196) 



For the potential 27ru(fc) = ^p- of a one dimensional 
plasma, if neglect the Landau damping and consider the 
long wavelength limit fc — >■ of the dispersion relation 
(|195p . we obtain at order k^: 



3(v^)k^ 



UJp 



Am? 



(197) 



This is the quantum generalization of the Langmuir wave 
dispersion relation. For the potential 27r{t(fc) = —2G/k'^ 
of a one dimensional gravitational system, we obtain 



= -2Gp + 3(w2)fc2 - (5(1.4) _ g^„2^2) ^ 



h^k 



2Gp Am^ 



(198) 



This gives the growth rate of the Jeans instability for small 
k when quantum effects are taken into account. Note that 
for a Fermi distribution at T = (see below), the term in 
parenthesis in equations (|197p and (|198p vanishes. 

Let us write the complex pulsation as uj = cOr + iiOi- 
For LUi — 0, the real and imaginary parts of the dielectric 
function e(fc, a;^) — er{k, Wr) + i^i{k, uir) are [66] : 



\k,uJr) = 1 - 2nu{k)P 



f{v)dv 



k I Am^ 



(199) 



27r^mu(A;) 

hk 



fik 
2m 



-f 



UJr 



hk 

2m 



(200) 



where P stands for the principal part. The imaginary part 
of the dielectric function is zero when = kvo where vq 
is solution of 



/(fo + hk/2m) = f{vQ - hk/2m). 



(201) 



If f{v) has a single maximum, then equation (|20ip has a 
unique solution vq. According to the Nyquist theorem, the 
system is stable ii er{k, kvo) > for all k and it is unstable 
if there exists some values of k for which er{k,kvo) < 0. 
Therefore, the mode k is stable if 



2TTu{k)P 



f{v)dv 



{v - vqY - 



>o, 



(202) 



and unstable otherwise. This generalizes the classical cri- 
terion [S3]. Equation (|202p can be rewritten [6 7) : 



1 + 2TTu{k) 



f{v„ + hk/2m) - f{v) 



{v - voY 



h?k^ 



dv > 0. (203) 



For the electrostatic interaction, and more generally for 
any repulsive potential u{k) > 0, a single humped distri- 
bution is always stable |67j . For a potential ti„ = — ^(Sn,±i 
corresponding to the (attractive) HMF model, a single 
humped distribution is stable iff 



f{v)dv 



l+'-P 

2 ./ {v~vor 



4 



> 0. 



(204) 



This generalizes the criterion obtained in |53| . 

As an illustration, let us consider the linear dynamical 
stability of the Fermi distribution at T = with respect 
to the Wigner equation. The DF is given by / = p/{2vf) 
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if < vp and f ^ it \v\ < vp. Substituting this dis- 
tribution function in equation (jl95p and solving for lj we 
obtain 



Am? 



9 9 Tik 

+ Vpk H Wfcotli 

m 



Tivpk 



2'Ku{k)pm 



.(205) 



For f{v) = pS{v), corresponding to vp — 0, this expression 
reduces to the form 



4™^ 



2TTu{k)pk'^ 



(206) 



It coincides with the fluid dispersion relation (jl75l) for a 
cold gas with Cs = 0. In the TF approximation ?i — >■ with 
Vp fixed, equation (j205p becomes 



„2 l2 



2TTu{k)pk'^ 



(207) 



It coincides with the fluid dispersion relation (|175|) with 
Cs = t;^ and ?i — > 0. Therefore, in the TF approximation, 
the dispersion relations of the Fermi distribution derived 
from the Euler and the Vlasov equations are the samej^. 
This is no more true when the quantum potential is taken 
into account (compare equations ()205p and (jl75p '). 

2 

For a potential 2iTu{k) = corresponding to a one 
dimensional plasma, we recover the results of [61j (the 
case of a gravitational plasma is obtained by reversing the 
sign of the interaction). On the other hand, for a poten- 
tial Un = — ^(5n,±i corresponding to the attractive HMF 
model (the repulsive HMF model is obtained by the sub- 
stitution k — > — fc), only the modes n = ±1 are allowed 
and the quantum dispersion relation (|195p becomes 



k 

1+2 



f{v)dv 



0. 



For the Fermi (or waterbag) distribution, we obtain 



(208) 



kM / 4TThVF\ 

17 \ kM J 



(209) 



where ${x) = x coth(.T) 
duces to the form 



For Vp ~ 0, this expression re- 



kM 
in 



(210) 



In the TF approximation ?),—>■ with vp fixed, we obtain 

(211) 



kM 
An ' 



This is the classical dispersion relation of the waterbag 
distribution It coincides with equation (|190p as ex- 
plained previously. This is no more true if we take into 



® Equation (|2U7p also represents the dispersion relation of 
the waterbag distribution in the classical limit. Therefore, in 
the classical limit, the dispersion relations of the waterbag dis- 
tribution derived from the Euler and the Vlasov equations are 
the same 1531. 



account the quantum pressure (compare equations (j209|) 
and (fT90)) '). For ?i ^ 0, using ^>(a:) = l + x'^ /3- x'^/45 + ... 
for .T — )- 0, equation (|209p can be expanded in powers of h 



4 



kM 
An 



n^fi^ 



4*8 

n n 



3fc2 45fc4 



.(212) 



This expression differs from equation (|190p by terms of 
order or smaller. In the limit h +oo, we get 



— + vp - hvp. 



(213) 



To appreciate the effect of the number of particles N, 
we introduce the dimensionless Planck constant ((77| and 
rewrite the pulsation equation (|209p in the form 



Anco^ 
kM 



Similarly, the pulsation relation (jl90p derived from the 
Euler equation can be written 



Anuj-^ 
kM 



2x! 

7V2 



1. 



(215) 



The thermodynamic limit corresponds to — > +oo with 
fixed X- In that case, the two relations take the form 



Anuj 
kM 



2^ 



1. 



(216) 



They return the critical Planck constant (|5T|). Therefore, 
the TF approximation is exact in the thermodynamic limit 
N +00. 

Remark: when coupled to a long-range potential of in- 
teraction, the classical Vlasov equation is known to de- 
velop filaments at smaller and smaller scales due to phase 
mixing and/or violent relaxation |36j . Considering (real 
or effective) quantum effects could be a way to regular- 
ize the Vlasov equation at small scales and serve as an 
alternative to coarse-graining. In that case, the (effective) 
Planck constant h could determine the scale at which the 
filaments are smoothed-out. 



6 Conclusion 

In this paper, we have generalized the HMF model in order 
to take into account quantum effects. We have considered 
the case of fermions and bosons with cosine interaction at 
T = 0. In the classical limit ?i — > 0, all the particles are 
located at = with velocity v = 0. This leads to a distri- 
bution function corresponding to a Dirac peak in position 
and velocity space: f{0,v) = M5{v)5{9). When quantum 
effects are taken into account, this Dirac distribution is 
regularized. Fermions tend to occupy the lowest energy 
states but cannot condensate into the state (0 = 0, = 0) 
because they must satisfy the Pauli exclusion principle. 
Bosons form a condensate in velocity space (all the bosons 
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are in the same quantum state with u = 0) but, because of 
the Heisenberg principle, they are delocalized in position 
space so that p{9) ^ M5{6). 

In the classical limit, the homogeneous phase is un- 
stable at T = 0. When quantum mechanics is taken into 
account, we find that the homogeneous phase can be stabi- 
lized when the normalized Planck constant x is sufficiently 
high. However, crucial differences exist between fermions 
and bosons. In the case of fermions, the stabilization is 
due to the Pauli exclusion principle while in the case of 
bosons, it is due to the Heisenberg principle. Therefore, 
the study of fermions is based on the Thomas-Fermi ap- 
proximation while the study of bosons is based on the 
Hartree approximation. This is why the thermodynamic 
limit N — >■ -f oo is different for fermions and bosons (see 
Appendix [X)). For fermions, the thermodynamic limit cor- 
responds to — > +0O in such a way that the coupling 
constant scale like k ^ N (yielding E ~ N^). This is 
a new prescription. For bosons, the thermodynamic limit 
corresponds to iV — > +oo with k ^ 1/N (yielding E ^ N). 
This corresponds to the Kac prescription like in the classi- 
cal regime. In the case of fermions, the stabilization of the 
homogeneous phase in the quantum regime occurs through 
a first order phase transition (as a function of the normal- 
ized Planck constant) while the phase transition is second 
order in the case of bosons. On the other hand, in the 
semi-classical limit ?i — ^ 0, the density profile of fermions 
is parabolic (see equation (|272[) ') while the density profile 
of bosons is Gaussian (see equation (|27ip ). 

One interest of the quantum HMF model is its rela- 
tive simplicity that should allow for accurate numerical 
simulations of quantum particles (fermions and bosons) 
with long-range interactions. In the present paper, we have 
only considered the thermodynamical equilibrium states 
of these systems, but the relaxation towards equilibrium 
is also of considerable interest. In particular, how does the 
system relax towards a steady state? What is the damping 
mechanism? These are important questions that we plan 
to address in the future. 

Acknowledgment: I am grateful to L. Dclfini for his 
assistance in some aspects of the numerical work. 



A Thermodynamic limits for fermions and 
bosons 

If we reintroduce the dimensional parameters, the Hamil- 
tonian of the HMF model reads 



1 ^ 

H = - mvl ~ gvn? ^ ( 



i=i 



R 



(217) 



where g is the coupling constant (the counterpart of the 
gravitational constant G in astrophysics), m is the mass 
of the particles and R is the system size (the particles 
are confined in the domain [— Tri?, 7ri?]). By comparing the 
different terms of the Hamiltonian, we get the scaling 



Therefore, the dimensionless energy is 



E 



gN'^m? 



(219) 



in agreement with equation ([25|). We can also define a 
dynamical time ^ R/v. Using equation (|218p . we get 



to 



R 



\J Ngm 



(220) 



Finally, in the classical regime, the temperature can be es- 



timated by ksT 



Using equation (j218p to estimate 



w, we obtain the dimensionless inverse temperature 

?y - PNgm^, (221) 

in agreement with the expression t] = /3fcM/(47r) of 

In the case of bosons at T = 0, the long-range interac- 
tion is balanced by the pressure arising from the Heisen- 
berg uncertainty principle AxAp ^ fi. Writing Ax ~ R 
and Ap ^ mv and using equation (j218p to estimate v, we 
obtain the dimensionless Planck constant 



XB 



^1/2^1/2^^3/2^ 



(222) 



in agreement with equation (|126l) . 

In the case of fermions at T = 0, the long-range inter- 
action is balanced by the pressure arising from the Pauli 
exclusion principle / ^ m? /h. Writing / ~ M/{Rvf) with 
vf ^ V and using equation (|218p to estimate u, we obtain 
the dimensionless Planck constant 



XF 



^gl/2^3/2' 



(223) 



E ^ Nmv'^ - N'^gm^. 



(218) 



in agreement with equation ([77)) . 

These scalings allow us to correctly define the thermo- 
dynamic limit. In the classical regime, the only dimension- 
less parameters are e and 77. The classical thermodynamic 
limit corresponds to iV ^ +00 in such a way that e and 
rj remain of order unity. If we take to ~ _R '--^ 1, we see 
that the choice g ~ 1/A^ (Kac prescription) leads to an 
extensive scaling E ^ N oi the energy and to an intensive 
scaling /3 ^ 1 of the inverse temperature. On the other 
hand, the dynamical time scales like tu ^ 1. We could 
also take g ^ 1 and m ^ This leads to the same 

scaling E/N ^ (3 ^ 1 oi the energy and inverse temper- 
ature but in that case Id ^ R/N^^^ (we may then take 
R ^ N^^'^ in order to have ~ 1). 

In the quantum regime, the dimensionless parameters 
arc e, 77 and x- The quantum thermodynamic limit corre- 
sponds to iV — >■ +00 in such a way that e, 77 and x remain 
of order unity. Let us take h ^ m ^ R ^ 1. In the case 
of bosons, we find from equations (|222[) that the coupling 
constant scales like g ^ 1/N which coincides with the Kac 
prescription. In that case, we have E/N ^ (3 ^ to ^ 1. 
Therefore, the bosonic quantum thermodynamic limit co- 
incides with the classical thermodynamic limit. In the case 
of fermions, we find from equation (|223p that the cou- 
pling constant scales like g ^ N which is different from 
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the Kac prescription. In that case, the energy scales hke 
E ~ which is different from the extensive scaHn^ 
and the inverse temperature scales like /3 ~ A^~^ which 
is different from the intensive scaling. On the other hand, 
the dynamical time scales like to ^ 1/N. Therefore, the 
fermionic quantum thermodynamic limit differs from the 
classical and bosonic thermodynamic limits. The fact that 
the bosonic and fermionic thermodynamic limits differ is 
not surprising since the quantum pressure which balances 
the long-range interaction corresponds to the Heisenberg 
uncertainty principle for bosons and to the Pauli exclusion 
principle for fermions which are of a completely different 
nature. 

We note that other scalings arc possible. For exam- 
ple, it makes sense to take ?i ~ m ~ 17 ~ 1 since these 
quantities are fixed physical constants that should not de- 
pend on A^. In the classical limit, this implies E ^ iV^, 
/3 ~ N-^ and to ^ R/Vn (we may then take R \/N 
in order to have ^ 1). In the case of fermions, using 
equations ((^^ and ^l^-l^^, we obtain R - iV^/^, 
E ~ 7V^, P ^ N^^ and to ~ 1. In the case of bosons, us- 
ing equations ([^^ and ([^ - ([2^ . we get R - N^'^/^, 
E N"^, (3 N-^ and to ^ N-^. We could also im- 
pose the Kac prescription R ^ m ^ gN ^ 1 to fermions 
by taking h N~'^. In that case E/N ^ j3 ^ to ^ I. 
Although this limit is mathematically conceivable, there 
is no physical reason why fi should depend on N . There- 
fore, we could impose the Kac prescription under the form 
h ^ m ^ gN ^ 1 by taking i? ^ A^. In that case 
E/N - ^ - 1 but - VN. 

These strange scalings should not cause surprise. They 
arise due to the long-range nature of the interactions mak- 
ing the energy non-additive. Similar (unusual) scalings are 
found for classical and quantum self-gravitating systems 
|l6l:j4l25| . 



We stress that these equations are equivalent to the mean 
field Schrodinger equation ((82|) -(|85 |) . In the classical limit 
?i — >■ 0, the quantum potential can be neglected and we 
recover the Euler equations at T = considered in [S]. 

Linearizing these equations around the homogeneous 
solution P = ^ I we obtain 



dSp 
'dt 



dSu 
'~d9 



dSu 



85$ 



6${9, t)^ / 5p{9\ t) cos{e - 9') d9'. 

271" Jo 



(227) 
(228) 

(229) 



Eliminating 5u between the first two equations, we get 



d^5p d'^S^ d^Sp 



89^ 



4 9614 



(230) 



Let us consider a perturbation of the density profile of the 
form 



Sp{e, t) = Re 



(231) 



Using equation (|229l) . the corresponding perturbation of 
the potential is 



, t) = Re 



kA 



(<5„,i+<5„,_i)e'("«-'^*^ 



(232) 



Substituting equations (j23ip and ()232D in equation (j230p . 
we obtain the dispersion relation 



LU = — n (d„.i + bn. 
47r 



(233) 



For n 7^ ±1, we get 



B Stability of the homogeneous phase for 
bosons 

In order to study the stability of the homogeneous 
phase, we shall use the Madelung representation of the 
Schrodinger equation in terms of hydrodynamical equa- 
tions. This will allow us to draw a close parallel with the 
stability analysis of the classical HMF model performed 
in [5]. For the bosonic HMF model, the quantum Euler 
equations ([III), ([Ml) and take the form 



(224) 
(225) 

(226) 



9p d , . ^ 
dt 09 ~ " 89 ~^ 2 d9\ ^ 



${9) = / p{9') cos{9 - 9') d9'. 

27r Jo 



The A'^^ scaling of the ground state energy can be directly 
seen on equation (|24|l . 



= — n^, (234) 

so that these modes are stable: the perturbation oscillates 
with a pulsation uj = fir? j^. For n = ±1, we obtain 



kM 
47r 



(235) 



In terms of the normalized Planck constant (|126p for 
bosons, we find that the homogeneous phase is stable iff 



X 



1/2 



> Xc 



72. 



(236) 



In particular, the homogeneous solution is unstable in the 
classical regime X — >■ [3] but it becomes stable in the 
quantum regime for x > Xc- In that case, it is stabilized 
by the Heisenberg principle or, equivalent ly, by the Bohm 
quantum potential. Note that the critical Planck constant 
Xc = 's/S precisely corresponds to the bifurcation point at 
which the inhomogeneous branch appears (see Appendix 
[C]) . For X < Xc, the inhomogeneous solutions with B ^ 
are stable while the homogeneous solutions are unstable. 
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C Bifurcation analysis for bosons 

In this Appendix, wc determine the behavior of the mag- 
netization and of the energy close to the critical point at 
which the bifurcation from the homogeneous (quantum) 
to the inhomogeneous (classical) phase takes place. 

The time independent mean field Schrodinger equation 
with a cosine potential can be written in dimensionless 
form (see Section H3|) : 



1 



X^^A" - b cos 0-ip = £ip, 



cos 9 dO, 



27r 



dO = 1, 



(237) 



(238) 



(239) 



(240) 



We assume that the wavefunction is symmetric with re- 
spect to = so that V'(~^) — Therefore, the 
boundary conditions are 



V;'(0) = i:\TT) = 0. 



(241) 



The homogeneous solution corresponds to V'o = 1/ \/27r 
and 6o = ^0 = 0. Close to the bifurcation, we make the 
expansion 



1 



/27r 



b = ebi+ (^b2 + £^63 + 



where e ^ 1 is a small parameter. 
At the order e, we obtain 



bi 



61 cos 6 = ]— £ii 



/2^ 



27r 



V'l cosO dd, 



(242) 

(243) 
(244) 
(245) 

(246) 
(247) 



At the order e^, we obtain 



b2 = 



2tt 



/27r 



cos 6 



V'2 + i't COS e d9, 



/27r 



2tt 



^2 de = — 



'2ti 



2ir 



i^l dO. 



(250) 
(251) 

(252) 



Solving these equations with the boundary conditions 
(PiT|) . we find that £2 = -&f/2, Xi = and 



V'2 



1 



169 cost 



/27r 8V27r 
At the order e'^, we obtain 



1 

&?cos(26l) (253) 



4V27r 



V-a - V2x2^i - (^biij2 + b2i>i + 



cos 9 



b.= 



27T 



£2lpl 



tp3 + 2V'ii/'2 COS 9 d9 



, (254) 



(255) 



Solving these equations with the boundary conditions 
(PiT]) . we find that £3 = -6162, X2 = -76?/(8V2) and 



36? 



16 



6162 



V2^V'3 = — + 63 cos 61 + — cos(26l) 



+^cos(30) + C, (256) 



where the constant C could be obtained by writing the 
normalization condition at order (but we shall not need 
it here). 

Combining the previous results, we find that the bi- 
furcation takes place at x = Xc = V^, i-e. at the point 
where the homogeneous phase becomes unstable. Close to 
the critical point, the magnetization, the energy and the 
wavefunction behave like 



7 



4\/2 



(Xc - x), 



(257) 



(258) 



2-K 



tpi d9 = 0. 



(248) 



Solving these equations with the boundary conditions 
((24T|) . wc find that £1 = 0, xo = V2 and 



■01 



/2n 



zbi cos 9. 



(249) 



[l + 6(x) cos 6*]. 



(259) 



We also find from our analysis that the total energy ctot 
vanishes at the order so that it scales like Ctot (xc — 
x)^. The prefactor could be obtained by extending the 
asymptotic analysis to next order. 
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D The semi-classical limit for bosons and 
fermions: harmonic approximation 

In the classical limit h ^ 0, the system is equivalent to a 
classical gas at T = 0. In that case, all the particles have 
collapsed at 6* = and the density profile is a Dirac distri- 
bution p{9) = MS (9) with magnetization 6 = 1. According 
to equations (|114p and ([7]), the total energy is etot = — 2 
and the energy appearing in the Schrodinger equation is 
£ = -1. 

In the semi-classical limit ?i — > 0, the distribution is 
strongly peaked around 9 — 0. We can therefore replace 
cos 9 by 1 — 9^/2 since 6* <C 1 and extend the integrals over 
the angles to infinity. This is valid in the limit /« — > 0. In 
that limit, equations (jl28p and (|129p become 



YV'" + ie'^ = (A + i)V, 



+ 00 



i;^ dO = 1. 



(260) 



(261) 



Using equation (jl39D . the total energy is 
etot ^-2 + 2x. 



(269) 



Finally, the wave function and the density profile are given 

by 



m = ( — 



1/4 



e 2x . 



and 



p{9) = TV — 



1/2 



(270) 



(271) 



These asymptotic expansions are compared with the nu- 
merical results in Section 14.51 

We can also consider the semi-classical limit for 
fermions. This corresponds to the limit x — >■ = 3/2 
in the equations of Section [3.4.21 For x = Xc, the density 
profile (|1D|) is a Dirac peak p{9) ~ MS{9). For x — > Xc, we 
can make the harmonic approximation cos 9 ~ 1 — 9^/2 
and we obtain 



For given k, equation ()260[) is just the ordinary 
Schrodinger equation of a harmonic oscillator. The fun- 
damental eigenvalue is 



X=-l + -K, 

and the corresponding wavefunction is 



1/4 

V'(0)=( — ) 

TTK , 



The magnetization can be approximated by 

r+oo / ra - 



d9 



(262) 



(263) 



(264) 



Substituting equation (|263p in equation (|264|) . we get 



V3M 



T\(x — Xc)^!"^ 



302 



4(x - Xc) 



1/2 



(272) 



for 9 <9c = V4(a;-a;c)/3 and p{9) = for 9c<9 <t^. 
In this limit, the integrals (j42p can be approximated by 



^'^"^''^ " 3^/2^'^ """^ 18V2 



m^{x-Xcf. (273) 



Therefore, when x Xc, the magnetization, the inverse 
polytropic temperature and the total energy are given by 



1 



[X - Xc), 
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2[x-XcY' 



(274) 



(275) 



b=l - -K. 
4 



(265) 



Finally, substituting the expansions (|262p and (|265p in 
equations (jl3ip and (|132p . we obtain 



X~K, £~-l+-K. 



(266) 



We can now express the results in terms of the normalized 
Planck constant x, defined by equation (I126L in the limit 
X — > 0. Eliminating k between equations (|265p and (|266p . 
we find that 



e ~ -2 + -{x - Xc). 



(276) 



Finally, using equations and ([7^ . the parameter x 
can be related to the maximal value of the distribution /x 
or to the Planck constant x by the relation 



1 



X 2{x-Xc)' 



(277) 
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